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I. Introduction 



The purpose of the present report is to describe in complete detail a FORTRAN code 
named Program SCAT 4 written by the UCLA group in order to analyze elastic scattering 
of various particles against complex nuclei by means of the diffuse surface optical model 
of the nucleus. 

While a number of similar programs have been prepared and used by other groups, 
there have been many requests for the UCLA program because of its flexibility and the 
availability of IBM 704 and 709 computers for which the program is written. 

The present program still contains some undesirable features and the UCLA group is 
constantly modifying it to make it more efficient and flexible. However, a "final" program 
will probably never be reached and it was decided to release Program SCAT 4 without 
further delay; as they develop, modifications and additions will be described in later 
reports. 

Other laboratories will probably add further modifications and the UCLA group will 
be grateful for description of such modifications as well as for any suggestions in this 
regard. Modifications and additions deemed worthwhile will be passed on to other users 
of the program but while the UCLA group is willing to serve partially as a central clearing 
house, the entire clerical responsibility cannot be assumed by the UCLA group. 

It should also be noted that, while every effort has been made to check out the program, 
the UCLA group cannot guarantee its complete correctness. 

Program SCAT 4 is available on a symbolic deck and will be mailed on request. Air 
mailing will require prepaid postage by requesting parties. 

Potential users of program SCAT 4 may find it useful to follow these suggestions in 
reading the present report: 

1) If the potential user is only interested in analyses with standard potentials he may 
proceed as follows: 

a) Read the introduction to the mathematical description. 

b) Consider the fundamental equations: (34), (35), (51), (78) through (85), (132), 
(137) through (139) in chapter II. 

c) Read chapter III, section A and the general flow chart. 

d) Read the description of subroutines INPT4 and OUTPT4 in chapter III, sec- 
tion B. 

e) Read chapter IV and VII. 

2) If the potential user is interested in all the features of the program, then a perusal of 
the whole report is advisable. The mathematical description of chapter II is a brief 
review of the theory and the basic equations are all listed there. Symbolic FORTRAN 
variables are indicated in capital letters and may be looked up in the glossary making 
up chapter V. 

Note that the program may be used for incident neutral particle by letting ZZ 1 = 0. 



II. Mathematical Description 



Program SCAT 4 calculates in the center-of-mass system the differential elastic scat- 
tering cross sections cr(9), the polarization P(0), and the total reaction cross section aji 
for particles of spin or 1/2 having any mass, charge and (non-relativistic) energy scat- 
tered by spinless nuclei of any mass and charge for various sets of diffuse surface optical 
model parameters. The incident and target particles are assumed to interact through a 
two-body potential consisting of a complex nuclear potential which includes spin-orbit 
interaction and whose shape can be specified by input parameters. When the incident 
particle is charged, the two body potential contains, in addition, the coulomb potential 
between an incident point charge and an extended, constant charge density target. 

The calculations include numerical integrations of the radial Schroedinger equations 
for the effective partial waves. The complex phase shifts are obtained as usual by matching 
the logarithmic derivatives of the numerically obtained nuclear wave functions to that of 
the coulomb (or spherical Bessel) functions. The phase shifts are then used to compute 
polarizations and cross sections which may be compared to the experimental values by 
means of the \ test. 

A. General Formulation 

We begin with a brief review of the basic theory relating to the scattering of spin 1/2 
particles by a zero spin target . We shall first consider the case of an uncharged incident 
particle and indicate later the modifications necessary if the incident particle is charged. 

The interaction is assumed to be of the form 

V T = V 1 + V 2 S-L (1) 

where V\ and V% are complex quantities depending only on the distance r between the 
incident particle and the target particle. In terms of the Pauli spin operator a, the spin 
operator of the incident particle, S, is given by 

S = \M (2) 

and the (relative) orbital angular momentum operator is given by 

L = rx(-f\. (3) 



The Schroedinger equation is then 

^ee J. Lepore, Phys. Rev. 79, 137 (1950). 



11 V 2 + V 1 (r) + V 2 (r)S-L 



^ = E^ (4) 



where 

V = : (5) 

mi + m b 

is the reduced mass, TOj and m^ being respectively the masses of the incident and target 
particles in atomic mass units. 

E = ~^E LAB (6) 

m i + m b 

is the energy in the center of mass system, -Blab being the lab energy of the incident 
particle in MeV. 

1. Uncharged Incident Particles 

The wave function corresponding to a wave incident in the positive z direction and 
normalized to one incident particle per unit time per unit area is 

* inc = -—e ikz Xinc (7) 



where v is the relative velocity, the wave number k is given by 



k = J 2 ^- = 0.2195376^/JiE ferrni" 1 (8) 

and the incident spin function is 

Xinc = a l/2 a + a -l/2^ ( 9 ) 

where a and (3 are normalized spin eigenfunctions of S z and a^ / 2 , o_i/2 the corresponding 
amplitudes. 

The partial wave expansion corresponding to (7) is given by: 



oc 



I Y^'~- --■'>■ " ' 47r 



*inc = ^E( 2£+1 )^( fcr )\/^^i y A^^ [°l/2« + «-l/2^J ( 10 ) 



where Ji(kr) is the regular spherical Bessel function of order £ and the normalized spherical 
harmonics are defined as 



m 



where Pi (cos 9) are the associated Legendre polynomials. 

The product functions YPa and Y\ r/3 which appear in (10) are simultaneous eigen- 
functions of the operators L , L z , S , and S z but not of the operator L ■ S which appears 
in the spin-orbit interaction. This may be remedied by introducing functions ty.* J which 



are simultaneous eigenfunctions of L , S , J , and J z and thus of L • S where J is the 
total angular momentum, 

J = L + S. (12) 

Since s = 1/2, the possible values of j are j = £ + 1/2 and j = £ — 1/2; the corresponding 
eigenfunctions are given by 



fW m 3 
^(.+1/2,1,. 



0)s rn i 

^e-i/2,e r 



2£+l 



' ' '"' ' ' /2 Y^- 1/2 a +J £ - m / + 1/2 YP +1/2 P, forj = l + l/2 



2£+l 



^-"^ + 1 /2 y m J -i/2 a + J / + ? + l/2 ^ +1/8ft for , = £ _ i/ 2 



2£ + 1 £ V 2£ + 1 

The incident wave function may now be written as 



(13) 



*i 




V 

4-7T 



J^V£+li e Je{kr) 




£=0 

oo 



a l/2 ^+1/2,^,1/2 + a -l/2 ^+l/2/,l/2. 



J] \/£ i e j e (kr) [-a 1 / 2 &t_ Wil 



&, 



-1/2 



1/2 l" a -l/2 ^-l/2,£,l/2 



(14) 



The total wave function can be written in a form similar to (14): 



* 



^ 



^«r 



total — Mnc "T" ^scatt 
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^ 



1/2 



^/ 



-1/2 




£=0 

oo 



1/2 ^+1/2^,1/2 ^ "-1/2 ^+l/2,€,l/2 



£ vs 



*7(»0 



.•£ * I 



1=0 



kr 



A/2 



1/2 



a l/2 ^-1/2,^,1/2 + a -l/2 ^£-1/2,^,1/ 



(15) 



where \l/j~ is the radial function associated with j = £ + 1/2 and ty „ is associated with 
j = £-1/2. 

The terms appearing in (15) are not coupled by the spin-orbit interaction, and sub- 
stitution into the Schroedinger equation (4) yields the following radial equations: 



dr 2 



2/i 



h 2 ( t 

V i + Y\ or ]l ~ 2 



^K = <> 



(16) 



where the quantity £ appears in the equation for ty » and — £ — 1 appears in the equation 
for $J. 

The radial wave function ^ must reduce to that of the incident wave, kr Ji(kr), when 
there is no interaction and must be such that only the outgoing wave is modified by the 
interaction. These conditions are satisfied by the asymptotic expression 



r± 



± 



Vf * krj e {kr) + Cf [-y e (kr) + ij e (jkr)) 



(17) 



which reduces to 



* 



± 



krj £ (kr) + C^e^ kr - £n ^ 



or equivalently 



^ ± = sin(fcr--) + C±e^-^/ 2 ) 



(" 



(19) 



as may be seen by applying the asymptotic expression for the regular and irregular spher- 
ical Bessel functions: 

kr jg(kr) = sin(fcr — £ir/ '2) 

kryi(kr) = — cos(kr — £ir/2). 



(20) 



c± 



On the other hand, in terms of complex phase shifts <L , (19) must be of the form 

$± ^ Af sin(fcr - £ir/2 + Sf) (21) 

Comparison of the coefficients of e and e~ in eqs. (21) and (19) yields 



Af = e i6 e 



2iS 



± 



1) 



(22) 
(23) 



Substituting (18) into (15) and subtracting ^ mc as given by (14), yields for ^ SC att the 
asymptotic form: 



*scatt 
where 



1 e 



ikr 



V r 



{A(0) 



a_l/ 2 e %<p a 



h/2 a + a -i/2l 3 \ +iB(9) 

B(B): 



a 1/2 e*P/3\ } (24) 



i 
k 



(25) 



£[C*-C7]P/(cos0) 
£=0 

The wave function of the scattered wave can more conveniently be expressed in terms 
of a and n, the unit vector normal to the scattering plane defined by 



ft sin 6 = k\ x k 







(26) 



where fcg an d k\ are unit vectors in the direction of propagation before and after scattering; 
thus 



^ s 



1 e 



ikr 



scatt 



V r 



[A(6) + B(9)a-n] X i 



1 e 



ikr 



V r 



-f(0)xi 



where f{6) is the operator 



f(6) = A(6) + B(6)o ■ ft. 



(27) 
(28) 



The differential elastic scattering cross section and polarization vector which are given 



by 



v(0) = (lf(0)XnJ[f(0)Xu 
P{9) = ± 



a{9) 



thus become 



a (6) = \A\ 2 + \B\ 2 + {A*B + AB*)n ■ P 



no) 



{\A 



2 -\B\ 2 )Pn 



A*B + AB* + 2\B\ 2 Pn-n 



n + i(A* B - AB*)n x P 



\A\ 2 + \B\ 2 + (A*B + AB*)P ■ n 
where the incident polarization vector Pq, is given by 

^0 = (4nc ^' 



(29) 
(30) 

(31) 
(32) 



(33) 



If the incident beam is unpolarized, i.e., Pq = 0, the scattered beam is polarized along 
the direction n, perpendicular to the scattering plane and 



a (9) = |v4| 2 + |£| 2 



(34) 
(35) 



Experimentally, the polarization is sometimes obtained from a double scattering ex- 
periment in the same plane wherein the polarization in the first scattering is known . 
The differential elastic scattering cross section for the second scattering may then be 
obtained from (31) and (35): 



\B\ 



1 



\A\ 



*2(6) = (\A\ 2 

= (\A\ 2 + \B\ 2 )(1 + P 2 -Pi) 
Referring to FIGURE 1, it is clear that 

m = n 2 



A*B + AB* _ 



\B\ 



n 2 -Pi 



(36) 



-nr 



•2. (37) 

so that the differential scattering cross sections along the r and £ beams are as follows: 

|2 , | D |2> 



a r 2 (0) = (\A\' + \B\')(l + P 2 P l ) 
a e 2 (9) = (\A\ 2 + \B\ 2 )(l-P 2 P l ), 



(38) 



2 L. Rosen, Proceedings of the International Conference on the Nuclear Optical Model, Florida State 
University, Tallahassee, 1959, pp. 72-90. 




Fig. 1 



the ratio of the scattering intensities becomes 



and solving for P2: 



which reduces when P\ = 1 to 






P2 



P2 



1 ^2- g 2 



a 2 -°2 



(39) 



(40) 



(41) 



2. Charged Incident Particles 

We next consider the case in which the incident particle has charge Ze and the target 
particle has charge Z'e. The potential V(r) must now include a term V c (r) which describes 
the coulomb interaction. For small values of r, V c will depend on the assumed charge 
distribution, while for large values of r, we must have 



ZZ'e 2 , , 
V c = (r large) 



(42) 



It is convenient to introduce the parameter r], 

7lJl 



V 



^ Z f e = 0,15805086 ZZ' 



&k ~~~ V #LAB 



m; 



(43) 



For the "incident wave" we take ^c( r )Xmci where ty c is the solution to the Schroedinger 
equation 



h ^0 ZZ'e 2 



2/x r 

corresponding to the scattering of two point charges. 



\fr r = E^ f 



(44) 



It is well known that in that case 



^ c = -^= T(l + i^e -1 / 2 ^ e iA;z F(-z?7, 1, ijfcf) (45) 



where £ = r — z and F is the confluent hypergeometric function. 

It is important to note that Sfr c includes a distorted incoming wave plus a scattered 
wave due to the point charge potential, and as such is not strictly an incident wave. 

The asymptotic form of \l/ c is given by 



*c 



_ I i[kz —r]£nk(r—z)} / i _ ^ 



where 



*V I V iA;(r-«) / (4g) 

— i (ft) p^kr— r\in2kr) 
r 



f(Q\ = 1 C -JV Msin 2 /2)+2i a /^x 



2ksm z ft/2 



is the Rutherford scattering amplitude and erg is given by equation (49), below, with 
1 = 0. 

The partial wave expansion of \I/ C is given by 



. e _ .^(a+D^S^/^W) (48) 

where ^(77, fcr) is the regular coulomb function and ag is the usual coulomb phase shift 
given by 

Of, = argT(£ + l + iT]) (49) 



Comparing equation (48) with (10) we see that in equation (14) it is necessary to 
^1- thus 

kr ' mUb ' 



replace ji(kr) by e tc7e 2' — -; thus, in this case 



oo 

■ 1/2 



£=0 



£ x/ZTl/e^ -*£-^ [ a 1/2 ^1/2,^1/2/2 + a_ 1/2 ^ +1/2>/jl/2 



(50) 




3^ 



— 2^ Vhe^ — |-a 1/2 ^_ 1/2A1/2 /2 + a_ 1/2 ^_ W)1/J 

£=0 

The total wave function can be written as a sum of the "incident" wave, \J/i nc , plus a 
"scattered" wave, \tscatt> where ^scatt now includes only interference terms and deviations 



from pure Rutherford scattering: 



* total = *i 



^scatt 




3^ 







£=0 



oc 



kr 



1/2 



1/2 %+l/2/,l/2l 2 + a -l/2 ^+1/2,^,1/ 



(51) 






£=0 



kr 



-0.1/2 %-l/2,£,ll^ 2 + a -V2 ^-l/2,£,l/2 



This wave function, ^ total > ls formally almost identical to the expression given by equa- 
tion (15) and the radial wave functions \&* obey an equation which is formally identical 
to equation (16) except that V\(r) must now include the coulomb potential V c (r) which 
may differ from a point charge potential at close distances. 

The radial wave function \& \ must now reduce to the "incident" wave, Fg{r), kr), when 
the potential becomes a coulomb point charge potential, and must be such that only the 
outgoing wave is modified by the non-coulomb interaction. These conditions are satisfied 
by the asymptotic expression: 



,± 



•y± 



*f * F £ (r), kr) + Cf [G £ (T], kr) + iF e (r], kr)] 



(52) 



which reduces to 

or equivalently 

^ = 



,± 



± i(kr-r)£n2kr-tir/2+a e ) 



sin 



*f 9±F e (r,,kr) + C?e 



(kr - V £n2kr - £ir/2 + a e ) + cfe i{ - kr ^ in2kr ' i ' K l 2+a ^ 



(53) 



(54) 



as may be seen by introducing the asymptotic expressions for the regular and irregular 
coulomb functions: 



Fi(r], kr) = sm(kr - r]in2kr - £ir/2 + c^) 
Gi(rj, kr) = cos(kr -rjin 2kr - £tt/2 + ag) 



(55) 



e± : 



In this case, the "nuclear phase shift" 5 P is taken to be such that the asymptotic form 



,± 



of \1/ p is given by 



4> 



± 



At sin(kr -rjin 2kr - £tt/2 + erg + St) 



(56) 



Comparison of the coefficients of e^ kr V? n2kr ) anc j e H fcr ntnkr) j n e q Ua tions (54) 
and (56) yields 



C 



± 



.4 



± 



1 r 

27 [ e 

J5± 



2i5i 



(57) 
(58) 
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Substituting (53) into (51) and making use of (46) and (50) we obtain for the asymptotic 
form of the total wave function 



* 



i i\kz—r) £nk(r—z)\ 
total = — ^= < e L 



1 



V 



Xi 



i A(kr—r\ £n2kr) . 
+ "-- " {^) 



y 



z&(r — 2;) 
l l/2 a + a -l/2^ +iB(6) 



—up 
a -l/2 e a ~~ a 



(59) 
« 1/2 e^] } 



where 



.. OO 



e=o 



B(6) = --J2e 2 ^[C+-Cl]P}(cose) 



(60) 



and /c(#) is given by equation (47). 

From this point, the formulation follows through as in the case of uncharged particles. 

B. Optical Model Potential 

1. Diffuse Surface Optical Model with Volume Absorption and 
Coulomb Spin-Orbit. 



The interaction (1) is assumed to have the form 

V T = VcN + ^SO + ^Coul + ^coul SO 



(61) 



where the terms appearing in equation (61) are respectively the central nuclear, spin-orbit 
nuclear, coulomb, and coulomb spin-orbit potentials. 

We shall first consider the case for which the real and imaginary parts of the central 
potential have a special common form factor (corresponding to volume absorption), and 
the spin-orbit potential is of the Thomas type. This particular central potential form 
factor has been used extensively and will be referred to as the standard form factor. We 
shall then discuss other form factors available in the program. 

(a) Central nuclear potential 



Vcs = (-V-iW)- 



(1 + e (r-R N )/aj \ ' 

where V and W are respectively the depths of the real and imaginary part of the nuclear 
potential in MeV (V and W are positive for an attractive, absorbing potential), and a 
common volume absorption form factor is assumed, where 



i?7V = RoN m h x 10 13 cm 



(63) 



i?ON being the nuclear radius constant and a is the rounding parameter in 10 cm 



13 



- 11 - 

(b) Nuclear spin-orbit potential 

The nuclear spin-orbit potential is often written in the Thomas form 



V S0 = X 



1 



1 d 



2M|c 2 1 r dr 



■V 



,(r-R N )/a 



1 



S-L 



(64) 



where M p is the proton test mass and c the velocity of light. If A were 1, the spin-orbit 
term would be that predicted by the Dirac equation. To provide more freedom in the 
model one writes 



A = 4 



m p\ V s + iW s 



M* 



V 



(65) 



where M^ is the pion rest mass and Vg and Wg are respectively the strengths of the real 
and imaginary parts of the nuclear spin-orbit potential in MeV. 

It may be noted that a negative value of the real part of A would be in accordance 
with the shell model of the nucleus where a (real) negative spin-orbit term is required to 
give the proper level sequence in contra-distinction to the atomic case. 



(c) Coulomb potential 

The coulomb potential is taken here to correspond to a constant charge density within 
the nucleus extending to a distance R c given by 



1/3 



R r = RnrmJ X 10" 1,3 CHI 



13 



where R nc is the coulomb radius constant; thus 



7'2 



2/ D 2>> 



Vcoul = (ZZ'e z /2Rc){3 - r z /R z c ) for r < R 



'2 



ZZ'e z /r 



for r > R ( 



(66) 



(67) 



(d) Coulomb spin-orbit potential 

The coulomb spin-orbit term is assumed to have the form 

1 d 



^Coul SO = (PP ~ 2> 



1 
Mp 2 



r dr 



^Coul 



S-L 



(6? 



where /dp is the proton magnetic moment in nuclear magnetons. It may be noted that 
the coulomb spin-orbit term is negligible except at very high energies. 

Substituting equations (62), (64), (67), and (68) into equation (16) and transforming 
to the dimensionless variable 

p = kr (69) 



3 W. Heckrotte, Phys. Rev. 101, 1406 (1956). 
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we find 



d 2 £(£ + 1) 
dp 2 p 2 

+ Lc) 


{ V T) 

rvs + iWs^ 

V £ ) 


( 


i ^ 




± e (p-PN)/kaJ 

Id/ 1 V 


pdp\i + e(P~PN)/ ka )_ 



e 

or 



±z 



-^Coul + ^CoulSO-U^(p) = (70) 



where 



^Coul 



Pc 
277/p 



3 - ^2 ) for P - Pc 



for p > p c 



U Covl SO 



1 f H \ 



2 V^"pc 
1/ h \ 2 



j {p P -\){2 n )(k 2 /p, 



fe) foP-$)M(*W 



or 

or 
-£-1 



and where 



Substituting now 



il^c 



P7V = fci?^ 
p c = kR c . 



2.00 x 10~ 26 cm 2 



2r] k 2 ■ \ 



h 



M P c 



2q 



E 



M P & 



2rj 



E 
931 



HP - \ = 2.7934 - 0.5 = 2.2934 
into equation (70) yields: 



p f{fi) 



-1 + 



£(£+1) 



V + iW 
E 



+ 



,{p~PN)/ka 



Vs + iW s \ (k\ (l 

E ) \ a ) \P (1 + e(fi~PN)/ka-j2 

£(£+1) (V + iW\ ( I 



1 + e (p-PN)/ka 



+ ^13 



for p < p c 
for p> p c 



0.004926 



r]E 
~PJ 



(71) 



(72) 



(73) 
(74) 



(75) 

(76) 

(77) 



or | > ^ff(p), for p< p c 



-1 + 



p 



E 



I _|_ e (p-p N )/ka 



+ 



2// 



+ 



y s + iW s \ ffc^ / 1 



{p-PN)/ka 



p (1 + e (p-Piv)/fea)2 



0.004926 



r/B 



or 
-£-1 



(78) 



*±(p), for p>p c 
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2. Nuclear Form Factors 

Equation (78) may be rewritten in such a way as to display explicitly the various 
nuclear form factors: 



d; 
dp 



2*?(p) 



11 



V S 2k f < \ _l ■ 



W s 2k 



V W 

p fcRKP) ~ ^ fci(p) • 
h hi p c 



fsi(p) - 0.004926^f 
PcJ 



I 

or 



11 I <?_ PI 

-2 
Pc 



,± 



tff (p), for p < p c 



i(£+l) V . . . .W , , , 2p. 

— 2 ^ /cr(p) " * ^ /ci(p) + - 



(79) 



Vq 2k , . , Wc 2A: . . «£ 

— - /SR P + *-# — /SI p) " 0.004926\ 

p° _ 



£■ a 



E a 



t 
or 



,± 



*f(p), for p > p c 



Three basic nuclear form factors and some special modifications of them are presently 
available in the program. In addition the coulomb spin-orbit term may be excluded at 
will. The required form factors may be chosen by assigning the proper values to the 
symbolic quantities KTRL as described on pages 33 ff. 

(a) Basic Form Factors 

(i) Volume absorption (ktrl(i) = o, i = l, 7, 8, 9, 10) 

1 



/cr(p) = /ci(p) 



/sr(p) = /si(p) 



(1 + e {p~PN)/ka^ 
1 e (p~PN)/ka 

p (1 + e {p-PN)/ka^2 



(80) 
(81) 



(ii) Gaussian absorption (ktrl(i) = l) 

/CR is g iven b y (80), /sr and / SI are given by (81) and 



fci(p) = e" 



w 



here 



-[{p-PG)/kbf 



1/3 



(82) 



PG = kRoG™<b", (83) 

Rqq being the nuclear Gaussian radius constant, and b determines the Gaussian width, 
(iii) Square well (ktrl(i) = 2) 

/CR(P) = /Cl(p) = 1 for P < PN 

= for p > p N 



Zsr(p) = /si(p) = 0. 



(84) 

(85) 
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(b) Special Central Nuclear Form Factors 4 

(KTRL(l) = o) 

The purpose of these form factors is to allow one to modify the knee or tail of the 
potential curve and produce central rises or depressions in the real and/or imaginary parts 
of the central nuclear potential, as specified by proper choice of the KTRL's. 

(i) FORM A (KTRL(7) = 1 for real part, KTRL(8) = 1 for imaginary part). 
/cr(p) and/or / CI (p) = [1 + h A (p)] f nAl {p) < p < p mA 



Pm A < P< PN 
PN < P < Pmax 



= fnAxip) 

= fnA 2 (p) 
(ii) FORM B (KTRL(7) = 2 for real part, KTRL(8) = 2 for imaginary part). 

/cr(p) and/or f CI (p) = [1 + h B {p)] f nBl (p) < p < p mB 

= fnB x (P) Pm B < P< PN 

= fnB 2 (p) PN < P < Pmax , 



36) 



(87) 



The presence of forms A and B allows distinct form factors in the real and imaginary 
parts. The presence of A\, A2 and B\ , B2 allows distinct shapes in the knee and tail of 
the form factors. Letting x be either A or B, and n be either nAi, 71A2, nB\, or nBy, 



h x (p) = 


= ho x 


k p ) 3 

\Pm x J 


\Pm x J 


+ 1 


= ho x ll- 


- P ) 2 ( 1+P ) 

Pm x J \ Pm x J 


(88) 


fn(p) = 


1 








(89) 


1+9n(p) 


where 












9n(p) = 


- exp < 


^n \ka ) 


XpnJ 


} 






(90) 



where Iiqa, Hq B , nA\, 71A2, nB\, nB^, Pm A , Pm B are selected constants. (The n's are 
always taken as > 0.) 

Note 1: If ho x is taken to be zero and nx\, nx2 are taken to be 1, forms A and B reduce 
to the volume absorption form. 

Note 2: The three curves defined by equations (86) and (87) join smoothly with con- 
tinuous derivatives as long as p mx is chosen less than p]\f. 

Note 3: Positive values of Hq x will produce central rises in the form factors while neg- 
ative values will produce a central depression. 



4 J.S. Nodvik, Proceedings of the International Conference on the Nuclear Model, Florida State Uni- 
versity, Tallahassee, 1959, pp. 16-23. 
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Note 4: If nx\ > 1, the knee of the potential will be sharper than for the usual volume 
absorption case, while < nx\ < 1 will soften the knee of the curve. 

Note 5: If n%2 > 1, this will shorten the potential tail while < 11x2 < 1 will extend 
it. 

Some typical shapes are presented in FIGURES 2, 3, and 4. 



(c) Special Nuclear Spin-Orbit Form Factors (ktrl(i) = 0) 

Two special nuclear spin-orbit form factors are available. They can be applied to the 
real and/or imaginary parts of the nuclear spin-orbit potential. The first of these form 
factors corresponds to the Thomas term applied to form A in the central nuclear potential, 
while the second uses form B itself; this permits one to study the result of deviations from 
the Thomas form. 

(i) DERIVATIVE FORM FACTOR A (KTRL(9) = 1 for real part, KTRL(IO) = 1 for imag- 
inary part) 



/sr(p) and/or f SI (p) = (ka) 
= (ka) 

= {ka) 
= {ka) 



1 d 
P dp 



(form factor A) 



1 dh A (p)\ (I df nAl 

-j fnA 1 (p) - (1 + h A (p)) — 

p dp J \p dp 

1 df nAl (p) 



(A 

dp 
for < p < p ma 



) (91) 



P dp 
1 df nA2 (p) 
P dp 



for p rrla < P< Pn 

for pN < P< Pmax 



where 



1 dh A (p) 6h 0A 
P dp 



Pm A 



Pm A 



1 df n p (p N \ 1 / p \ n 2 

p dp V ka/ p z \pn J 



(92) 
(93) 



and fn(p) and gn{p) are given by equations (89) and (90). 

(ii) FORM FACTOR B (KTRL(9) = 2 for real part, KTRL(IO) = 2 for imaginary part) 

1 



/srG°) and/or /gi(p) = -■ [form factor B as per equation (87)] 



(94) 



Note: If Jiqa ls taken to be zero while nA\ and nA2 are taken to be 1, the derivative 
form factor in (91) becomes identical to the usual spin-orbit form factor (81). 
Some typical shapes are presented in FIGURES 5, 6, and 7. 
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3. Final Formulation for Machine Calculation 

The complex radial wave function ^ f (p) may be written as 



Vf(p) = xf(p) + iyf(p) 



(95) 



and equation (79) for a ■ £ = £ or — £ — 1 can now be separated into two real coupled 
differential equations, and dropping the subscripts and superscripts for convenience: 



d z x 

= px - qy 

d 2 y 

— = qx+py 



(96) 



where 



P = ^CR + ^SR ( or 

-£-1 



£{£+!)) 



q = U C l + U S i or 



(97) 



Formulas (97) are convenient for programming purposes as the f/'s are now independent 
of I. indeed: 



^CR 



!- p/CR+— (3-^0 ) forp<p c 

E Pc\ pi 



for p > p c 



= 


-1^? 


U C I = 


- ¥ /ci 


U S R = 


Vq 2k riE 

---§-— fSR- 0.004926-^ 

E a p% 




Vq 2k riE 

= -#— /SR- 0.004926^ 

E a p 6 




bi E a bi 



(98) 



(99) 
(100) 
(101) 



4. Numerical Integration 

Equations (96) must be integrated numerically twice for each £ = to £ m ax where 
£ max +l corresponds to a partial wave negligibly disturbed by the scattering. 

The method chosen for numerical integration is the 3-point Runge-Kutta method: it 
lends itself to easy starting, permits one to change the interval quite easily and gives 
excellent accuracy with relatively large steps. 



20 



fc 






2 
U 



CO 






Q- 



Hh 



is 

fa 






o 

O 



> 



■3 
o 

T3 



3 

-: 

CO 



o 

o 



" <J < < < 

< « n n n 

is is £ £ 

^ b b b b 







Pi 

.2 

+j 

.3 

'in 

X 

> 



$ 



o 
c 

T3 



o 

.3 



r& 



CO 



o 



ISjIjI io/pire hSjM 



21 



fc 



LO 

CO 






O 

o 
.03 



o 

'in 

o 

T3 



a 
Pi 

CO 



o 



< 

fa 



fa 






lO - 









fa 



CO 

LO 



O 

o 

< 

fa 









i — 1 


e@ 


©(§ 


fa 







1 


1 


j 


i 


i 


1 


- 


1 








- 


- 




^A 


q\ 




- 


- 


©/ 


e 




@A 


- 


/ 1 


1 


i > 


i 


i 


i\ 






o 



oo 
o 



OJ 

o 



o 
o 



o 
o 

I 



OJ 

o 
o 

I 



oo 
o 

o 

I 






3 

,o 






id 
=1 

fa 



LO 



.a 



3 

fa 

I 

O 

id 






3 
O 



LO 



fa 
(H 
O 

fl 

'fa 

CO 

I 

CO 



fa 



ISJJ lo/pire hSjM 



22 



fc 









CO 









Cm 



CM 
Z 

fa 



IN 

in 

o 



o 






o 



O 
T3 



3 

-: 

CO 



£3 

O 

o o o 



www 



o 
© 



<;<;<; 



fa 



fa 



13 

fa 









- 




- 


- 




- 


@> 


co) 


V 


- 








- 









LO 



-* 



Q. 



3 

o 



3 



3 
> 

C 

o 

Pi 



CO' 



a 
> 

'in 

o 

o 

+j 
o 
.a 



.o 



OJ 



t-l 

o 

I 

a 
"3, 

CO 



o 



CO 

o 



o 



ISjIjI io/pire hSjM 



23 



Given xn, y a , in, m, at p h where i n = \J*) etc. 



xn = f{xa,yn,pi); m = g{xn,yn,Pi) (102) 

Ap Ap 

Xi2 = Xii+xn — ; Vi2 = yn + m-^- I 103 ) 

Xi2 = f(xi2,Vi2,P%+-^-)\ yi2=9(Xi2,Vi2,Pi + ^) (104) 

(Ap) 2 (Ap) 2 . . 

Xi3 = Xi2 + xn - ; 2/^3 = Vi2 + j/ii (105) 

Ap Ap 

^3 = /W3,Z/i3,Pi + ^-); j/i3 = #W3,Z/i3,P* + -^-) (106) 

• Ap .. (Ap) 2 . Ap .. (Ap) 2 ,_ 

zi4 = z»2 + :Eti-2- + Zi2— g— ; yj4 = yi2 + %i^- + ?/i2^— (107) 

^4 = f(xn,yiA,Pi + Ap); y i 4 = g(x H ,y i 4,p i + Ap) (108) 

and finally 

arj+1,1 = Xji + Axi = x& -\ — {xn + x i2 + ^3) + Ap % (109) 



Xi+1,1 = in + Aii = in + —(in + 2£ i2 + 2x i3 + i lA ) (110) 

J/t+1,1 = V%\ + Ayj = 2/ii + (j/ji + 2/i2 + y»3) + Ap^i (111) 

2/i+l,l = 2/ti + A^ = j/ii + -g-(y*l + 2i/i2 + 2^3 + ^4) (112) 

The process is continued until the nuclear potential becomes negligible at which time 
the wave functions and their first derivatives must be saved for later matching with those 
of the coulomb function. 

Starting values: If Pj n iti a i is very small, the following starting values may be used: 

Xt(p = Pinitial) = (Apl) m ; X t (p = pinitial) = (^ + l)(Api) . , 

Ve(p = Pinitial) = 0; Mp = Pinitial) =0 

5. Coulomb Functions 

The regular and irregular coulomb functions are given by the following asymptotic 
formulas which may be used successfully for large values of p: 



Fo~sin[Re(v?o)]e" Im( ^ o) 1 
Fi^sintReOpOle- 1111 ^ 1 ) 
G ~cos[Re(^ )]e" Im(w) 
Gi-costRe^i)^- 1 ™^ 1 ) 



(114) 



where 
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ipo = p — r) £n 






(115) 



and where 



a\ = -n, a2 



a k 




(116) 



with a similar recurrence formula holding for b^ 



a = argT(l + in) 
<J\ = uq + tan - n 



(117) 



Furthermore the quantity ao may be successfully approximated over the whole range of n 
by the following formula: 



^0 



-n 



7 -l en 



\ ln(77 2 + 16) + -tan _1 ( 



12(n 2 + 16) 



1 n 2 -41 



-1 -1 fl 

tan ?7 + tan I - 



1 n A - 160?7 2 + 1280 



tan" 



-l en 



(T) 



(118) 



30 (?7 2 + 16) 2 105 (I6 + 77 



2\4 



The above formulas which can of course be generalized for any value of £ are equivalent 
though not formally identical to the formulas listed by Abramowitz and by Froberg . 

Rather than use these formulas for obtaining Fg and Gg for any value of £ > 1, it is 
preferable to make use of recurrence formulas. 

The following upward recurrence formula is suitable for finding Gf 



G 



(2£ + 1; 



£+1 



1(1+1) 



G £ -(£ + l)[£ 2 +n 2 ] 1/2 G e 



£[(£ + !] 



n 



21 1/2 



(119) 



Tables of Coulomb Wave Functions, Vol. I, National Bureau of Standards, Applied Mathematics 
Series 17, Washington, 1952, p. XV. 

6 C. E. Froberg, Rev. Mod. Phys. 27, 399 (1955). 
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A similar recurrence relation can only be used for downward recurrence on the Fgs, 
otherwise results rapidly lose all significance. This may be done by means of a method 
due to Stegun and Abramowitz and which is essentially as follows. 

Let it be required to compute Fg from £ = to £ = £ majX - 

(1) Let tO) = £ max + 10 

(The number 10 is arbitrary but has found satisfactory from practical experience) 



,(1) 



,(1) 



,(1) 



Let F\A = and F\A =0.1. Successive values of F« can be computed from 
£ = to £ = (S > — 1 by means of the downward recurrence formula: 



F, 



(1) 



(2£ + i; 



i(i+i) 



F, 



(1) -i[(i+l) 2 + V 2 ] 1/2 F^ 



£-1 



(£+l)[£ 2 + v 2 ] 1/2 
Letting the constant 

a = (F^ ) G 1 -F^G )(l+r ] 2 ) 1 / 2 
one may compute successively 

Ft 
for £ = £ m ^ + 1 to £ = 0. 



(120) 



Fi%-> 



(121; 



(122) 



(2) To verify the accuracy of the Fg's obtained above one may compute as above a new 
set of functions F\ starting perhaps from £^ > = V* > + 5 (again the number 5 is 



obtained from practical experience) and letting now F 



(2) 



yields a new set of F£s. 



£( 2 ) + l 



0, F 



(2) 



^(2) 



0.1. This 



(3) Comparison of the two sets of F^s obtained in (1) and (2) above indicates the 

accuracy of the computation. If this proves insufficient, let £^ > = £^ > + 5 and 

(3) (3) 

starting from F ,1 = 0, ^/ 3 ^ =0.1 one may obtain a third set set of Fg's which 

is to be compared with the second set. 

This procedure may be continued until two successive sets of Fg 7 s are found to agree. 
The derivatives of the coulomb functions may be obtained from the formula 



(*+l) 



V 



Y e -[(£ + l) 2 + V 2 ] L ^Y ( 



21 1/2 



£+1 



Y, 



where Yd stands for either Fo or Go. 



V 



(123) 



e 



( 



7 Stegun and Abramowitz, Phys. Rev. 98, 1851 (1955). 
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6. Phase Shifts 

The phase shifts are obtained in the usual fashion by matching the logarithmic deriva- 
tives of the coulomb functions with those of the numerically integrated functions at a 
value of p sufficiently large so that the nuclear potential becomes negligible. 

Matching the logarithmic derivative of the nuclear function \H^ = X£ + iy£ with that 
of its asymptotic form 

F i + (G e + tF e )C e 



yields 



which lead to 



* e = F e + (G e + iF e )C e 
^ F e + (G e + iF i )C e 



C 



± 



*P'i 



*?F t 



*fG l 



*P't 



i(¥± F e - Vffy 



(124) 



(125) 



the quantities Cg being related to the complex phase shifts through equation (57). 

7. Cross Section and Polarization 

The differential elastic scattering cross section a(8) and the polarization P{8) for an 
unpolarized incident beam are obtained from equations (34) and (35) while the reaction 
cross section may be obtained as follows. 



°R 



N ahs 



(126) 



where iV a k g is the absorbed flux, and N- mc is the incident flux which was assumed to be 1 
(see equation (7)). By definition, 



n. 



abs 



h 
2i\x 



* 



t 3^ total 



total 



Ol 



* 



d^ 



total " 



t 
total 



Ol 



n 

r sin 8 d8 dip 



(127) 



where the integral is taken over the surface of a large sphere of radius r = tq. Substituting 
equation (51) for ^ total ^o equation (127) and making use of the orthonormality of the 



W- / 's and of the relation 



l l/2 



a -l/2 



1- 



yields after carrying out the surface integration: 



-*-»*.-$±w{>>{-£) 



*P d 
kr dr 




4tt 

V 






h 

2ip 



V d (*\ 



kr dr \ kr 



*l d 



V 



kr dr \ kr 



kr dr \ kr 



(128) 



r=r 

(129) 



r=r 
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Now substituting the asymptotic form (52) for ^! . and making use of the Wronskian 
relations 

G £ F £ ~ F £ G £ = 1 
we are led to the following: 



47T 

V 



h 

2ifi 



^* d /*£ 



*? d /*^* 



kr dr V kr ) kr dr \ kr 



4tt 



r=r 



¥ [lm(C±) - |C±| 2 ] 



(130) 
(131) 



Finally, substitution of (131) into (129) yields 



4n 



£=0 



° = vi E W + x ) Im ^ + ) - (mw - ( Re (^ + ))' 



Im(C-) - (im(Cr)) 2 - (Re(C7)) : 



(132) 



Note: The quantities e l appearing in equation (60) may be obtained by the following 
recurrence formulas: 

~ (i + l) 2 -V 2 } 
— -s o cos 2a > 

_(£ + l) 2 +r] 2 

\£ + l) 2 -ri 2 



Re(e 2iai + 1 ) = cos2a £+1 

Im(e 2 ^+ 1 ) = sin2a m 



sin 2a i 



_(( + i) 2 + v 2 

while the Legendre polynomials obey the usual relations 



2ry(l + l) . 

-s ^- sm 2a* 

(£ + l) 2 +?7 2 £ 

2?](£ + 1 



L(£ + l)2 +r? 2 



COS 2(7/ 



(133) 



P (cos 0) = 1 , Pi (cos 0) = cos 9 
P m (cos0) = TTT [(2£+i)coB9Pf(coB0)-ePi- 1 (coa9)] 



pP(cos 



£ + 1 



sin 6? 



cos #P^(cos 6) - P e+1 (cos 6)] 



One may also compute the Rutherford scattering cross section: 

a c (6) = \f c (6)\ 2 . 



(134) 
(135) 

(136) 



8. Chi Square Deviation 

Experimental and theoretical quantities may be compared by means of the chi square 
deviation: 

XT = X 2 a + Xp (137) 

where 

" a th (#)-a ex (#) n ' 



e e 



Aa ex (9) 



>th 



p u \e) - P ex (e) 

AP ex (9) 



(138) 
(139) 



-28- 

where the a (9) and P (9) are the theoretically obtained cross sections and polarizations 
while a ex (9), Aa ex (9), P ex (#), AP ex (9) are respectively the experimentally given cross 
sections, standard deviations in the cross sections, polarization and standard deviations 
in the polarization. 

It should be noted that the constants were chosen such that the differential and reaction 
cross section will be obtained in units of 10 cm . The polarizations are of course 
dimensionless ratios. 

9. Normalization 

The radial wave functions ^ -, and their derivatives obtained from numerical integra- 
tion of the radial Schroedinger equation contain an arbitrary normalization factor, l/M„ . 
This factor however does not affect the cross section and polarization since these are ob- 
tained from the phase shifts which in turn are obtained from ratios of logarithmic deriva- 
tives (see equation (125)) wherein the Mg's cancel out. If on the other hand the normalized 
radial wave functions and their derivatives are required, the normalization terms may be 
obtained as follows: 

The asymptotic form of ty \ must obey equation (52) but improper normalization 
results in the fact that the calculated wave functions are actually given by 

xf(p) + iyf(p) = Mf {F e ( v , p) + Cf [G £ (rj, p) + iF e ( v , p)] } (140) 

Now, for p < p max the nuclear potentials are negligible and equation (52) represents the 
exact solution; in particular, at p = p m ax, we must have 

x £ (pmax) + iyfiPm&x) = Mf {F e (i],p max ) + Cf [G £ (r], p max ) + iFifa pmax)]} (141) 



whereby 



= _^ 4(j>m*x)+ivtU..~,_ _ {[[2] 



± _ X £ (pmax) + iy £ (pr 

Fe(V> Pmax) + C £ [G e {T], p max ) + iF e (r), pmax)] 



and the normalized radial wave functions and their derivatives are given by 

-ip- = Mj [x z {p)+iy z {p)] 



(143) 



and the complete normalized wave function is given in equation (51) with vJ ' » as above 
in equation (143). 

Note: During the numerical integration the program may renormalize the wave func- 
tions and their derivatives at any value of p for which overflow takes place by dividing the 
functions and their derivatives by the largest of these. This is accompanied by an explicit 
printout as explained in the description of subroutine RKINT. Such occasional internal 
renormalization must of course be taken into account if correctly normalized functions are 
required. 



III. Program Description 



A. General Description 

1. Machine Specifications 

Program SCAT 4 has been written for an IBM 704 with floating point traps or an 
IBM 709, with a 32,768 words memory, no drum and a minimum of two tape units. 

The program can probably be modified for a 16K memory by reducing the number 
of 0's (up to 75 allowed here) and the number of ts (up to 50 allowed here). A large 
part of the memory (7500 words) is occupied by the Legendre polynomials and this may 
also be reduced by computing the polynomials as required. Furthermore, the program 
contains a large number of printouts which may be abbreviated to save storage space. 

2. General Program Description 

The program was designed to compute cross sections, polarizations and chi square 
deviations at a number of specified points in the space of the optical model parameters V, 
W, A, VS, WS, and if needed BG (RO, RC and RG are kept fixed), for a given set of input 
data. 

The time to carry out a run for a single set of parameters depends of course upon 
the maximum values of £ and p; for p-Cu at 10 MeV (£ m ax = 10, p = 0.0625 (.0625) 0.50 
(0.25) 10.0) a run takes about 45 seconds including about 15 seconds for maximum output 
to tape. 

The program has been written in the form of subroutines to allow easy checking and 
modification. Some of these subroutines are not yet available, but some provision have 
been made to include them in the future. The following subroutines written in FORTRAN 
are specific (sub)routines of the program: 



Main routine 


MAIN4 






Subroutine 


CTRL4 






Subroutine 


INPT4 


Subroutine 


PGEN4 


Subroutine 


POT1CH 


Subroutine 


INTCTR 


Subroutine 


POP1 


Subroutine 


RKINT 


Subroutine 


SIGZRO 


Subroutine 


CSUBL 


Subroutine 


FSUBC 


Subroutine 


AB 


Subroutine 


EXSGML 


Subroutine 


SGSGCP 


Subroutine 


RHOTB 


Subroutine 


SIGMAR 


Subroutine 


COULFN 


Subroutine 


CHISQ 


Subroutine 


RMXINC 


Subroutine 


OUTPT4 



The following subroutines are general utility routines used by the program: 
Subroutine - SKIP written in FORTRAN 
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Subroutine - LEAVE written in FORTRAN 

Subroutine - SPILL written in FAP 

The following subroutines are used in conjunction with the Load-and-Go system in use 
at WDPC (Western Data Processing Center, UCLA). The effect of using this system is 
described in section III-A-3 below. 

Subroutine - SAVE 

Subroutine - PDUMP 

Subroutine - EXIT 



The program assumes the presence of the following Fortran elementary function subrou- 
tines: 



LOGF 


(natural logarithm) 


SINF 


(sine) 


COSF 


(cosine) 


EXPF 


(exponential) 



SQRTF (square root) 

ATANF - (arc tangent) 

3. Use of the WDPC Load-and-Go System 

Program SCAT 4 has been written for the Load-and-Go system in use at the WDPC, 
UCLA. This only affects it as follows: 

(i) Special subroutines of the load-and-go system. 

Subroutine SAVE 

The purpose of this subroutine is to allow the operator to interrupt the calculation 
without loss. The program is normally run with Sense Switch 1 off; turning on Sense 
Switch 1 will cause the program to call SAVE after completing the innermost DO loop 
of subroutine CTRL4. SAVE then writes on tape the content of the core memory as well 
as all other information required to continue the computation such as the contents of the 
AC, MQ, index registers, etc. . . . 

A restart routine will then later reload the core memory, reset all registers etc. . . , and 
return right after the CALL SAVE statement. The following statements up to statement 
number 66 are then required to properly position the input data tape as the latter was 
probably rewound when the computation was interrupted. 

To eliminate the use of subroutine SAVE, remove from subroutine CTRL4 all statements 
from statement number 118 to statement number 66 inclusive. 
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Subroutine PBUMP(a,(3) 

The purpose of this subroutine is to provide a partial core dump of all quantities 
between the location of the arguments in the call statement. Subroutine PDUMP is called 
by subroutine LEAVE whenever difficulties such as overflow or division by zero take place. 

To eliminate subroutine PDUMP, replace in subroutine LEAVE the statement CALL 
PDUMP(A,ZZ) by whatever statements will cause the required core dump. 

Subroutine EXIT 

This subroutine terminates the program. 

To eliminate subroutine EXIT, replace statement number 151 in subroutine INPT4 by 
whatever statement will be used to terminate the program. 

(ii) end Statements. 

The usual FORTRAN END statements do not appear in the program as the load-and-go 
system provides them automatically. 

(hi) Input and Output Statements. 

In conjunction with the load-and-go system, the program is input from tape, while 
the input data is brought in from tape 7 and all the output is to tape 6. 

All these particular features can of course be easily modified to use the program either 
directly or in conjunction with any other system. 

4. Error Indications: 

(i) Division by zero. 

Every division which could conceivably have a zero divisor either because of the range 
of numbers used or because of an error in the input data is followed by an IF DIVIDE CHECK. 
Detection of a zero denominator is then followed by an explicit print out and a CALL LEAVE 
statement which leads to the next set of input data. In order to be sure that no division 
by zero remains undetected, every subroutine which contains an IF DIVIDE CHECK state- 
ment also begins with an IF DIVIDE CHECK to verify that the trigger is off at the start 
of the subroutine; if the divide check trigger is found on at the start, there is an explicit 
printout to that effect followed by a CALL LEAVE statement. 

(ii) Overflow. Underflow. 

Overflow and underflow are monitored by subroutine SPILL (JSPILL, ISPILL, x, y) which 
needs only be called once by MAIN4. When SPILL is called, it replaces the quantities JSPILL 
and ISPILL by zeros. Thereafter, in case of overflow (underflow) the subroutine replaces 
the overflowed (underflowed) quantity with x (y) and places into JSPILL (ISPILL) the 
address of the command which caused overflow (underflow) to occur for the first time. 
Program SCAT 4 uses x = y = 0. 

Every subroutine in which computations are carried out starts by setting ISPILL and 
JSPILL equal to zero to insure correct identification of possible subsequent overflow or 
underflow. The subroutine then ends with a check of ISPILL and JSPILL. If either of these 
is not zero, there is an explicit printout describing the overflow or underflow. Underflow 
results therefore in substituting zero for the underflowed quantity, but the computation 
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proceeds. Overflow on the other hand results in substituting zero for the overflowed 
quantity and leads to a CALL LEAVE statement to stop the computation. 

B. Detailed Descriptions of the Specific Routines of the 
Program 

MAIN4 

The main routine which is only used at the start of the program carries out the 
following steps: 

1) Calls SPILL which controls overflow and underflow (see III-A-4-ii). One such call 
statement is sufficient to put SPILL in permanent control for all subroutines. 

2) Sets up EPS1, EPS2, EPS3, which are constants used to control the accuracy of the 
Coulomb functions computations, and EPS4 which is used in subroutine POT1CH. 

3) Inputs identification and program numbers. 

4) Calls CTRL4. 

CTRL4 (Control 4) 

This subroutine controls the whole flow of the program. It was coded as a subroutine 
to allow it to be called by subroutine LEAVE. It carries out the following steps: 

1) Advances group identification and resets run identification numbers. 

2) Call INPT4. 

3) Calls POT1CH. 

4) If KTRL(5) = 1, calls POP1 
if KTRL(5) = 0, proceeds. 

5) Calls SIGZRO, FSUBC, EXSGML. 

6) Sets up five (or six) nested DO loops for varying successively V, W, a, V s , W s (and 
b for a surface absorption potential). The following steps are always done within 
the innermost DO loop: 

a) If Sense Switch 1 is on, calls SAVE 
if Sense Switch 1 is off, proceeds. 

b) Advances run identification number. 

c) Calls RHOTB, COULFN, RMXINC, PGEN4, INTCTR, CSUBL, AB, SGSGCP, SIGMAR. 

d) If KTRL(2) = 0, proceeds 

if KTRL(2) = 1, calls CHISQ. 

e) Calls OUTPT4. 

7) When all the DO loops have been completed, returns to step 1. 
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INPT4 (Input 4) 

1) Inputs KTRL(l); if KTRL(l) = 100, calls EXIT 

if KTRL(l) ^ 100, proceeds. 

2) Inputs KTRL(I), I = 2 to 13. 

3) Inputs FMI, FMB, ELAB, ZZ, RC, V, W, RO, A, VS, WS, RG, BG, DV, DW, DA, DVS, 
DWS, DBG, HA, PMA, FN1A, FN2A, HB, PMB, FN1B, FN2B, NVMAX, NWMAX, NAMAX, 
NVSMAX, NWSMAX, NBGMAX. 

4) Sets up TV = V to TBG = BG (starting values of the parameters). 

5) Inputs NMAX, forms NMAXP = NMAX-1. 

6) Inputs RHOIN(I), I = 1 to NMAX and DRHOIN(I), I = 1 to NMAXP. 

7) Computes FMU as per equation (5) 
Computes ECM as per equation (6) 
Computes FKAY as per equation (8) 
Computes RHOBN as per equation (73) 

Computes RMA and RMB (see Glossary, under PMA, PMB) 
Computes RHOBC as per equation (74) 
Computes ETA as per equation (43). 

8) Inputs LMAXM, forms IMAX = LMAXM + 1. 

9) Sets IIN(J) = 1, J = 1 to LMAX (see description of subroutine INTCTR) 

10) If KTRL(5) = 0, proceeds 

if KTRL(5) ^ 0:a) inputs JMAX 

b) inputs THETAD(I), I = 1 to JMAX 

c) computes THETA(I), I = 1 to JMAX. 

11) If KTRL(2) = and/or KTRL(3) = 0, proceeds, 
if KTRL(2) ^ and KTRL(3) ^ 0, inputs 

SGMARX(I), DSGMEX(I), POLEX(I), DPOLEX(I), I = 1 to JMAX. 

12) Returns to CTRL4. 

POT1CH (potential 1 check) 

The purpose of this subroutine is to check whether £ max is sufficiently large so that all 
the partial waves sensibly affected by the potential are included and to check whether p ma x 
(the point at which the coulomb functions will be matched to the nuclear wave functions) 
is sufficiently large to insure that the non-coulomb part of the potential is negligible. If 
^max and/or p max are too small, the subroutine increases them, and sets IIN(£ max )= 1- 
The quantities p max and £ max may be checked or not according to the value assigned to 
KTRL(13): 
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KTRL(13) = 1: check both £ max and p m ax 
KTRL(13) = 2: check p max only 
KTRL(13) = 3: check £ max only 
KTRL(13) = 4: do not check either. 

Pmax and £ m ax are checked in various ways depending upon the potential form. The 
routine operates as follows: 

1) The routine first calculates the maximum values of V, W, A, VS, WS, and, in the 
case of a Gaussian absorption, of BG over the specified grid of these parameters. 

2) If KTRL(l) = 0, standard potential (or variation thereof), the routine checks, if re- 
quired, that: 



a) Pmax is sufficiently large so that 

(^2 + ^2)1/2 



< e 4 . (144) 



E (1 + e (Pmax-P7v)/fca) 

If this condition is not met, p max is increased by the last value of Ap and the 
check is repeated. This is accompanied by the print out: 

RHOIN(NMAX) = (value of old p max ) + (last value of DRHOIN) 

RHOIN(NMAX) IS TOO SMALL IN NUCLEAR POTENTIAL, 
b) The routine also checks, if required, that £ m ax is sufficiently large so that 

W 2 + W 2 1 



E (1 + e (4nax-pAr)/fca) 



< e 4 . (145) 



If this condition is not met, £ max is increased by 1 and the check is repeated; 
this is accompanied by the following printout: 

LMAXM = (value of old LMAXM) + 1 

LMAXM TOO SMALL BECAUSE OF CENTRAL POTENTIAL. 
The routine then checks that £ m ax is sufficiently large so that 



., v% + w% i 

2k 2 ^y_^ & i < (146) 

E (1 + e (4nax-/07v)/fca) ~ V ' 

If this condition is not met, £ m ax is increased by 1 and the check is repeated; 
this is accompanied by the following printout: 

LMAXM = (value of old LMAXM) + 1 

LMAXM TOO SMALL BECAUSE OF SPIN ORBIT POTENTIAL. 

3) If KTRL(l) = 1, Gaussian absorption, 
a) The check on p max is as follows: 

< e 4 ; (147) 



E (1 + e (Pmax-PAr)/fca) 
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and 

E e -(Pmax-PG/^) 2 < e4 . (14 8 ) 

E 

If these conditions are not met p max is increased as before and the checks are 
repeated; this is accompanied by the same printout as above. 

b) The check on £ m ax is as follows: 

VI , , 

Tf 7^ wi < e 4! (149) 

and 

^ e -(W-p G A&) 2 < e4 (150 ) 

and as in equation (146). 

If these conditions are not met £ m ax is increased by 1 and the checks repeated. 
The prints-out are given on the previous page. 

4) If KTRL(l) = 2, Square well 

a) The check on p max is as follows 

Pmax > PN (151) 

b) The check on £ m ax is as follows 

W > PN + 3. (152) 

Failure to meet these conditions leads to increases in p m ax and/or £ max accompanied 
by the same printouts as given above, after which the checks are repeated. 

The program uses EPS4 = 0.001. This quantity is specified in the MAIN4 routine. 
The checks described above are based on a rough estimate of the phase shifts using a 
WKB expression. 

POP1 

Computes P(L,J), PP(L,J), L = 1 to LMAXP, J = 1 to JMAX as per equations (134) and 
(135) and returns to CTRL4. 

SIGZRO (Sigma zero) 

Computes SIGMAO and SIGMA1 as per equations (117) and (118) and returns to CTRL4. 

FSUBC 

Computes FCR(J) and FCI(J), J = 1 to JMAX as per equation (47) and returns to CTRL4. 

EXSGML (Exponential sigma £) 

Computes EXSGMR(J), EXSGMI(J) for J = 1 to LMAX as per equation (133) and returns 
to CTRL4. 
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RHOTB (Rho tabulation) 

The purpose of this subroutine is to construct a table of p's and Ap's corresponding to 
each step of the numerical integration. This table is formed from the arrays of RHOIN(I) 
and DRHOIN(I) which are input by subroutine INPT4 



Input 


Arrays 


RHOIN(I) 


DRHOIN(I) 


RHOIN(l) 


DRHOIN(l) 


RHOIN(2) 


DRHOIN(2) 


RHOIN(NMAX-l) 


DRHOIN(NMAX-l) 


RHOIN(NMAX) 





Computed Tables 


RHO (I) 


DRHO(I) 


RHO(l) 


DRHO(l) 


RHO(2) 


DRHO(2) 


RHO(ILAST-l) 


DRHO(ILAST-l) 


RHO(ILAST) 





p= RHOIN(l) (DRHOIN(l)) RHOIN(2) . . . (DRHOIN(NMAX-l)) RHOIN(NMAX) 

RHO(I+l) = RHO (I) + DRHO(I) 

DRHO(l) = DRHO(2) = ■ ■ ■ = DRHO(I) = DRHOIN(l) 

up to RHO(I) = RHOIN(2), etc 

RHO(l) = RHOIN(l); RHO(ILAST) = RHO(NMAX) 
LLAST>NMAX. 

If RHOIN(NMAX) is given in such a way that it cannot be reached by an integral 
number of DRHO(I)'s, the last interval is shortened (up to 50%) or lengthened (by no 
more than 50%) so that RHO(ILAST) = RHOIN(NMAX). 

COULFN (Coulomb functions) 

This is the most complex subroutine of the program. It computes the regular and 
irregular coulomb functions and their derivatives for L = 1 to LMAXM at p= RHOMAX by 
means of asymptotic formulas. The main steps are as follows: 

1) The a and b series appearing in equation (115) are calculated according to equa- 
tions (116) and are cut off when either: 

(a) The term N a (or Nfr) is such that the next term exceeds in magnitude the 
previous one, i.e., when 

[Re(U Na + l)] 2 + [lm(U Na + l)] 2 > [Re(U Na )} 2 + [Im^J] 2 (153) 



where 



U k 



a k 



(k 



i\ Jfc-1 

-UPmax 



(154) 



and similarly for the b series. 

(b) The contributions of both the real and imaginary terms give undetectable con- 
tributions to the real and imaginary parts of <pq (and similarly for ipi). During 
these computations, the value of p m ax may be increased by addition of the last 
value of DRHOIN and the computation starts all over again under the following 
condition: 
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a) The a or b series is identically equal to zero. This is accompanied by the 
printout: 

SERIES IN PHIO OR PHI1 IS ZERO, CHECK DATA, IF OK 

INCREASE RHOMAX = (value of old RHOMAX) + (value of last DRHOIN) 

b) Either of the two series diverges too quickly, i.e., the iV a -th (or A^-th) term 
still gives a non-negligible contribution to the series obtained so far, viz. 



[Mu Na )\ 


2 + 


Mu Na )] 2 




iRe^t^Uk)] 


2 
+ 


^(eis 1 ^); 


2 



> EPS3 



(155) 



(EPS3 is given the value 0.00001 in the MAIN4 routine.) 
This is accompanied by the printout: 

IF OK A OR B SERIES DIVERGES TOO QUICKLY 

INCREASE RHOMAX = (value of old RHOMAX) + (value of last DRHOIN). 

c) Over 48 terms are required in either the a or b series. This is accompanied 
by the printout: 
INCREASE RHOMAX = (value of old RHOMAX) + (value of last DRHOIN) 
A OR B SERIES CONVERGES TOO SLOWLY. 

2) The quantities (po, ipi, Fq, F\, Gq, G\ are formed according to equations (114) and 
(115), and the Wronskian is checked for accuracy requiring that 



i + v z 



-1/2 



Fod - F X G Q 



i + v z 



-1/2 



< EPSl 



(156) 



(EPSl is given the value 0.00001 in the MAIN4 routine.) 

If this condition is violated p max is increased and the computation starts all over 
again; this is accompanied by the following printout: 

INCREASE RHOMAX = (old value of RHOMAX) + (last value of DRHOIN) 

BAD INITIAL WRONSKIAN. 

3) The regular coulomb functions are formed by downward recurrence as per equa- 
tions (120) and (122) according to the accompanying description. 

Agreement between successive sets of Ffs is verified by checking that 



(An) /jp( n + 1 ) 



1 



< EPS2 



(157) 



(EPS2 is given the value 0.00001 in the MAIN4 routine) for £ = to £ max . 

During this computation the value of p max is increased and the computation starts 
all over if it turns out that £(i) > 4nax + 40. This is accompanied by the printout: 

INCREASE RHOMAX = (old value of RHOMAX) + (last value of DRHOIN) 
L TOO LARGE IN FBAR(L). 
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4) The irregular coulomb functions are formed by upward recurrence as per equa- 
tion (119) and the Wronskian for every £ = to £ m ax + 1 is checked for accuracy 
requiring that 



F e G e+i - F i+i G e 



£ 



[(£ + l) 2 + V 2 ] 



21 V2 



< EPS1 (15? 



(EPSl is given the value 0.00001 in the MAIN4 routine.) 

If this condition is violated the value of p m ax is increased and the computation starts 

all over again; this is accompanied by the printout: 
INCREASE RHOMAX = (old value of RHOMAX) + (last value of DRHOIN) 
BAD WRONSKIAN FOR L = (value of £ + 1 for which equation (158) failed). 

5) Finally the derivatives of the coulomb functions for £ = to £ max are formed as per 
equation (123). 

RMXINC (Rho max increase) 

The purpose of this subroutine is to extend the table of RHO(I) and DRHO(I) by in- 
crements of the last value DRHOIN until the final value of RHO(I) equals RHOMAX which 
may have been increased by the subroutine COULFN. 

PGEN4 (Potential generator 4) 

The purpose of this subroutine is to form tables of the ^-independent parts of the 
potential corresponding to the RHO(I) tables and suitable for using in the numerical inte- 
grations. 

These include: 

UCRB(I), UCIB(I), USRB(I), USIB(I) for I = 1 to ILAST and corresponding to the values 
at the beginning of an interval of integration; a corresponding table of form factors is also 
formed: 

FFCR(I), FFCI(I), FFSR(I), FFSI(I), 
and 

UCRM(I), UCIM(I), USRM(I), USIM(I), 
and 

FFCRM(I), FFCIM(I), FFSRM(I), FFSIM(I) for I = 1 to ILAST - 1 corresponding to the 
values in the middle of an interval of integration. 

The original and tightest part of the subroutine corresponds to a standard form factor; 
modifications have been added to permit use of a variety of form factors briefly described 
earlier. 

The subroutine operates as follows: The UCR-'s are calculated as per equation (98), 
the UCI-'s as per equation (99), the USR-'s as per equation (100) and the USI-'s as per 
equation (101), wherein: 

(i) KTRL(I) = 0: VOLUME ABSORPTION OR SPECIAL NUCLEAR FORM FACTOR: 
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If KTRL(7) = 0, 

= 1, 
= 2, 



If KTRL(8) = 0, 

= 1, 
= 2, 

If KTRL(9) = 0, 

= 1, 
= 2, 

If KTRL(IO) = 0, 

= 1, 
= 2, 

(ii) KTRL(l) 



/cr is computed as per equation (80); 
/cr is computed as per equation (86); 
/cr is computed as per equation (87); 

f(jl is computed as per equation (80); 
f(jl is computed as per equation (86); 
f(jl is computed as per equation (87); 

/sr is computed as per equation (81); 
/sr is computed as per equation (91); 
/sr is computed as per equation (94); 

/si is computed as per equation (81); 
/sr is computed as per equation (91); 
/sr is computed as per equation (94); 

i: Gaussian absorption 

/cr is computed as per equation (80); 
f(jl is computed as per equation (82); 
/sr is computed as per equation (81); 
/si is computed as per equation (81); 

(iii) KTRL(l) = 2: SQUARE WELL 

/cr is computed as per equation (84); 
f(jl is computed as per equation (84); 

/sr and /ci are taken to be zero. 
Furthermore, 

If KTRL(ll) = 1, USR- are computed as per equation (100) including the coulomb 
spin-orbit term. 

If KTRL(ll) = 0, USR- are computed as per equation (100) excluding the coulomb spin- 
orbit term, i.e, the second term on the right hand side. KTRL(7) to KTRL(ll) can of course 
be given any combination of permitted values. 

INTCTR (Integration Control) 

For each value of L = 1 to LMAX this subroutine carries out the following steps: 

1) Sets up starting values for the numerical integration as per equation (113). The 
quantities IIN(L) are not especially useful at the present time, but they have been 
included in order to permit start of the numerical integration at various values of p 
depending on £ and thus permitting considerable time saving by foreshortening the 
numerical integrations. A study of this method is presently under way. 

2) Calls RKINT which performs the numerical integration. 



[FFCR] 8 


= /CR 


[FFCR] 


= /CR 


[FFCR] 


= /CR 


[FFCI] 


= fci 


[FFCI] 


= fci 


[FFCI] 


= fci 


[FFSR] 


= fsR 


[FFSR] 


= /srA° 


[FFSR] 


= /sr/2 


[FFSI] 


= /si 


[FFSI] 


= /siA° 


[FFSI] 


= /Si/2 


[FFCR] : 


= /CR 


[FFCI] : 


= fci 


[FFSR] : 


= fsR 


[FFSI] : 


= /si 


[FFCR] : 


= /CR 


[FFCI] : 


= /ci 



8 FFCR refers to the symbolic variables FFCR(I) and FFCRM(I) appearing in the program (see glos- 
sary of symbols), similarly for FFCI, FFSR, and FFSI. 
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3) Stores the final values of the functions and their derivatives at the completion of 
each integration. 

RKINT (Runge-Kutta integration) 

This is the most crucial subroutine in the program as most of the time is spent in nu- 
merical integration. Special efforts have therefore been made to produce a rapid program. 

The subroutine integrates numerically as per equations (102) to (112) the differential 
equations (96) operating simultaneously on the two sets corresponding to a ■ £ = £ and 
-£-1. 

Special provisions have been made to avoid overflow; this is accomplished by dividing 
all the functions and their derivatives by the largest of these at every step (RENORM); 
whenever such renormalization is carried out it is accompanied by the following printout: 
RENORMALIZATION FACTOR = (value of RENORM) IN RKINT FOR CODED 
L = (value of £ + 1) and RHO = (value of p at which renormalization took place). 

CSUBL 

This subroutine computes Cp as per equation (125) for £ = to £ ma , x - 

AB 

This subroutine computes A(J) and B(J) for J = 1 to JMAX i.e., for the various angles 9 
required, as per equation (60). 

SQSGCP (Sigma, sigma-coulomb, polarization) 

This subroutine computes a(6>), P(6), <J C {9), as per equations (34), (35); (136) and 
finally a(6)/a c (6) for the various angles required. 

SIGMAR 

This subroutine computes (tr as per equation (132). 

CHISQ (Chi Square) 

This subroutine computes Xa(@)i Xai Xp(#)? Xp? Xj> as P er equations (137), (138) 
and (139). 

Note: The quantities Aa ex (6) and AP ex (0) are always assumed to be non-zero. Thus to 
avoid including an unknown experimental quantity, the corresponding standard deviation 
must be taken as very large. 

OUTPT4 (Output 4) 

Several output formats are available: 

(1) Minimum output (KTRL(6) = l). 

(a) Basic quantities 
NUMPRG 

KTRL(I) for I = 1 to 13 

FMI, FMB, ELAB, ZZ, V, W, A, RO, VS, WS, RC, BG, RG RHOBN, RHOBC, 
RHOBNG, ECM,ETA, FKAY, FKAYA, FKAYB 
and, if either KTRL(7), (8), (9), or (10) is not zero, 
HA, RMA, FN1A, FN2A, PNA, HB, RMB, FN IB, FN2B, PMB, 
then RHOMAX, LMAXM, NMAX, RHOIN(I) for I = 1 to NMAX, 
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DRHOIN(I) for I = 1 to NMAX-1, SGMRTH 
and, if KTRL(2) = 1, CHI2ST, CHI2PT, CHI2T. 

(b) Basic Table 

THETAD(I), SGMATH(I), SRATIO(I), POLTH(I), 

and, if KTRL(2) = 1, SGMAEX(I), POLEX(I), for I = 1 to JMAX. 

(2) Normal output (KTRL(6) = 0) 

(a) Basic quantities 
(See above) 

(b) Basic Table 
(See above) 

(c) Form factor table (output only if KTRL(12) = l) 
RHO(I), FFCR(I), FFCI(I), FFSR(I), FFSI(I), 

for I = 1 to ILAST. 

(d) Fitting table (output only if KTRL(2)=l) 
THETAD(I), DSGMEX(I), DPOLEX(I), CHI2S(I), 
CHI2P(I), CHI2(I) for I = 1 to JMAX. 

(e) L table 

L, CR1(L), CI1(L), CR2(L), CI2(L) for L = 1 to LMAXM (corresponding to £ = 
to tniaxj- 

This output is made for every run, and maybe preceded by underflow descriptions 
which may be ignored, and by other comments referring to an increase in p m ax, ^max, 
renormalization, etc. 

Every page of output is headed by the run number on the left and the page number 
on the right. The number of lines per page is held to be less than 50, otherwise the 
subroutine calls subroutine SKIP which starts a new page. 

SKIP 

This subroutine increases the page number, resets K, the line counter, and outputs 
the run and page number. Note that arguments giving the number of lines, page and run 
numbers are required. 

LEAVE 

This subroutine is called whenever a run gets into difficulty because overflow, or divi- 
sion by zero occur. The subroutine calls PDUMP to give a partial core dump. 

This subroutine was included so as to allow for various possible requirements upon 
overflow and division by zero without having to change every command where the difficulty 
might occur. 



IV. Description of Input Data 



All data is input from tape 7. The input data tape is prepared from IBM cards which 
contain one piece of input data per card in either of the two following formats: 



Columns. 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


1415 


72 


Integers 
Floating nos. 


X 

± 


X 




X 


X 
X 


X 
X 






















X 


X 


X 


X 


X 


X 




± 











fractional part 



Note: Any floating point format which uses 15 columns or less and is acceptable to 
FORTRAN may be used in place of the above. 



(1) The following identification data is input first: 



NUMRUN(l) 
NUMRUN(2) 
NUMRUN(3) 
NUMRUN(4) 
NUMRUN(5) 
NUMPRG 



month 

day 

year 

set number (put in to start with 1) 

run number (put in to start with 1) 

program number (we use 4). 



Note: The identification which consists of the five quantities NUMRUN(I), I = 1 to 5, 
is printed at the top left of every output sheet. NUMRUN(4) is advanced every time 
a new set of data is input, NUMRUN(5) is advanced every time a run is made with a 
new set of parameters. 

(2) Then, for every set of run, i.e., for every set of input data: 



(a) Controls 

KTRL(l) = 
= 1 
= 2 

KTRL(2) = 
= 1 



Standard potential (possibly with generalized form factors) 
Gaussian absorption 
Square well 
no x required 
X required 



9 The quantity A is eventually discarded but it must still be input as 1/2 to avoid overflow in the 
early part of the program. 
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KTRL(3) = 


= 


same experimental values as in last set 




= 1 


new experimental values coming 


KTRL(4) 




not used in present program 


KTRL(5) = 


= 


same angles as in last set 




= 1 


new angles coming 


KTRL(6) = 


= 


normal output 




= 1 


minimum output 


KTRL(7) = 


= 


UCR - Standard form 




= 1 


UCR - form A 




= 2 


UCR - form B 


KTRL(8) = 


= 


UCI - Standard form 




= 1 


UCI - form A 




= 2 


UCI - form B 


KTRL(9) = 


= 


USR - derivative standard form 




= 1 


USR - derivative form A 




= 2 


USR - form B 


KTRL(10) = 


= 


USI - derivative standard form 




= 1 


USI - derivative form A 




= 2 


USI - form B 


KTRL(11) = 


= 


do not include coulomb spin-orbit 




= 1 


do include coulomb spin-orbit 


KTRL(12) = 


= 


do not print out form factors 




= 1 


do print out form factors 


KTRL(13) = 


= 1 


check p max and £ max 




= 2 


check p max only 




= 3 


check £ max only 




= 4 


do not check p max nor £ max . 



(b) Basic data 

FMI, FMB, ELAB, ZZ , RC, V, W, RO, A, VS , WS, RG, BG, DV, DW, 
DA, DVS, DWS, DBG, HA, PMA, FN1A, FN2A, HB, PMB, FN1B, FN2B, 

NVMAX, NWMAX, NAMAX, NVSMAX, NWSMAX, NBGMAX. 

(c) Integration data 

NMAX, RHOIN(I) for I = 1 to NMAX, DRHOIN(I) for I = 1 to NMAX - 1, 

(d) LMAXM 

(e) Angles: 

if KTRL(5) = 1 input: JMAX, THETAD(I) for I = 1 to JMAX 

(f) Experimental data: 
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KTRL(3) = 1 also requires KTRL(2) = 1 for proper operation. 
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if KTRL(2) = 1 and KTRL(3) = 1 input: 



SGMAEX(I) 

DSGMEX(I) 

POLEX(I) 

DPOLEX(I) 

(3) Final card: 
KTRL(l) = 100. 



for I = 1 to JMAX 
for I = 1 to JMAX 
for I = 1 to JMAX 
for I = 1 to JMAX 



V. Glossary and Description of Symbolic 
Variables Appearing in Common and Dimension 

Statements 



FORTRAN Symbol 


Math. Symbol 


Description 


A 


a 


Rounding parameter appearing in stan- 
dard potential, see eq. (62) 


AR(I), AI(I) 
I = 1 to 75 


Re{aj}, Im{aj} 


1) Real and imaginary parts of the 
terms of the auxiliary series used to 
calculate asymptotically the coulomb 
functions, see eq. (116) 




Re{A(9i)}, Im{A(^)} 


2) See eq. (60) for definition 


BR(I), BI(I) 
I = 1 to 75 


RejfeJ, Im{6j 
Re{B(0i)}, lm{B($i)} 


1) Ibid, see eq. (116) 

2) See eq. (60) for definition 


BG 


b 


Width parameter in Gaussian absorp- 
tion see eq. (82) 


CHI2(I) 
I = 1 to 75 


x 2 m 


= xl(0i) + x 2 P (0i) 


CHMP(I) 

I = 1 to 75 


xpm 


See eq. (139) 


CHI2PT 


x 2 P 


See eq. (139) 


CHI2S(I) 
I = 1 to 75 


x 2 M) 


See eq. (138) 


CHI2ST 


Xa 


See eq. (138) 


CHI2T 

CR1(L), CI1(L) 
for L = 1 to 51 


x 2 

Re(C+), Im(C+) 


= xl + Xp 

See eqs. (57) and (125) 


CR2(L), CI2(L) 


Re(Cj-), Im(Gy) 


See eqs. (57) and (125) 


DA, DV, DW, 
DVS, DWS, DBG 




Amount by which A, V, W, VS, WS, 
BG must be incremented for succeed- 
ing runs (these increments may be in- 
put as positive, zero or negative). 


DPOLEX(I) 
for I = 1 to 75 


AP ex (^) 


Standard deviation in the experimental 
polarization (must never be input as 0) 


DRHO(I) 

for I = 1 to 250 


Ap* 


Interval of numerical integration (see 
description of subroutine RHOTB) 


DRHOL 




Last interval to be used in the numeri- 
cal integration 
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FORTRAN Symbol 


Math 


. Symbol 


Description 


DRHOIN(I) 
I = 1 to 250 






Interval of numerical integration spec- 
ified by input for RHOIN(I) < p < 
RHOIN(I+l) (See description of subrou- 
tine RHOTB) 


DSGMEX(I) 
I = 1 to 75 


Aa ex (^) 




Standard deviation in the experimen- 
tal differential elastic scattering cross 
section in square fermis/sterad, (must 
never be input as 0) 


ECM 


E 




Incident energy in center-of-mass sys- 
tem (MeV) 


ELAB 


-^LAB 




Incident energy in laboratory system 
(MeV) 


EPS1, EPS2, EPS3 


ei, £2, £3 




Error thresholds appearing in various 
parts of the calculation of the coulomb 
functions. See eqs. (155) to (158) 


EPS4 


H 




Error threshold used in POT1CH sub- 
routine, see eqs. (144) to (150) 


ETA 


V 




See eq. (43) 


ETA2 


2 






EXSGMR(L), 
EXSGMI(L) 
L = 1 to 51 


Re{e 2 ^}, Im{e 2ia e} 


See eq. (133) 


F(L), L = 1 to 52 


Ft 




Seeeq. (114) and (122) 


FBAR(L), 
L = 1 to 91 


F (n) 




See eq. (120) 


FCR(I), FCI(I) 
I = 1 to 75 


Reifcie,)}, lm{f c (9i)} 


See eq. (47) 


FFCR(I), 
FFCRM(I) 
I = 1 to 250 


fCR(pi) 
fcR.(Pi + 


2 ) 


Form factors for the real central part 
of the potential at the beginning and 
middle of an integration interval (See 
eqs. (80), (84), (86), (87) and descrip- 
tion of subroutine PGEN4) 


FFCI(I), FFCIM(I) 
I = 1 to 250 


fci(pi) 

fCR(Pi + 


&Pi\ 
2 ) 


As above for the imaginary central part 
of the potential (See eqs. (80), (82), 
(84), (86), (87), and description of sub- 
routine PGEN4) 


FFSR(I), 
FFSRM(I) 
I = 1 to 250 


Isr(Pi) Isr(pi + ^r) 


As above for the real spin-orbit part 
of the potential (See eqs. (81), (85), 
(91), (94) and description of subroutine 
PGEN4) 
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FORTRAN Symbol 




Math. Symbol 


Description 


FFSI(I), FFSIM(I) 


fsi(Pi) fsiiPi + tP) 


As above for the imaginary spin-orbit 


I = 1 to 250 






part of the potential (See eqs. (81), 
(85), (91), (94), and description of sub- 
routine PGEN4) 


FKAY 


k 




See eq. (8) (inverse fermis) 


FKAYA 


ka 






FKAYB 


kb 






FMB 


m b 




Mass number of target nucleus (atomic 
units) 


FMI 


mi 




Mass number of incident particle 
(atomic units) 


FMU 


H 




Reduced mass of incident particle 
(atomic units (see eq. (5)) 


FN1A, FN2A 


nA\ 


nA 2 


See eq. (86) and following description 


FN1B, FN2B 


nB\ 


nB 2 


See eq. (87) and following description 


FF(L), L = 1 to 51 


*i 




See eq. (123) 


G(L), L = 1 to 52 


Gg 




Seeeq. (114) and (119) 


GP(L), L = 1 to 51 


G e 




See eq. (123) 


HA, HB 


hoAi 


h QB 


See eq. (88) 


IDATA 






Number of sets of data to be processed 
after making use of subroutine SAVE 


IFIRST 






Initial value of I, the subscript appear- 
ing in RHO(I) 


ILAST 






Final value of I, the subscript appear- 
ing in RHO(I) 


IIN(L), L = 1 to 51 






Originally designed to allow input of 
any desired value of IFIRST for various 
L's in order to speed up the numerical 
integration. In the present program the 
IIN(L) are all set equal to 1 by subrou- 
tine INPT4 


ISPILL, JSPILL 






Underflow and overflow indicators used 
in conjunction with subroutine SPILL 


JMAX 






Total number of angles input 
(JMAX<75) 


JMAXT 






Temporary storage for JMAX used after 
calling subroutine SAVE 
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FORTRAN Symbol 


Math. Symbol 


Description 


KTRL(I) 




Controls used throughout the program 


I = 1 to 13 




to specify the potential, input and out- 
put type (see description of input data) 


KTRLT(I) 




Temporary storage for KTRL(I) used af- 


I = 1 to 13 




ter calling subroutine SAVE 


L 


£ + 1 




LMAX 


tmax + J- 




LMAXM 


tmax 




NA, NV, NW, 




DO loop variables used in subroutine 


NVS, NWS, NBG 




CTRL4 to specify the number of times 
the parameters have been incremented 


NAMAX, NVMAX, 




Total number of incrementations of the 


NWMAX, 




parameters specified as input data ( > 


NVSMAX, 




1) 


NWSMAX, 






NBGMAX 






NINPUT 




DO loop variable used after calling sub- 
routine SAVE in order to count the 
number of sets of processed input data 


NMAX 




Total number of input values of 
RHOIN(I) specified in input 


NMAXT 




Temporary storage for NMAX used after 
calling subroutine SAVE 


NMAXP 




= NMAX - 1 


NUMPRG 




Program number (see description of in- 
put data) 


NUMRUN(I) 




Identification (see description of input 


I = 1 to 5 




data) 


POLEX(I) 


P ex (0 t ) 


Experimental value of the polarization 


I = 1 to 75 






POLTH(I) 


P th (0i) 


Calculated value of the polarization See 


I = 1 to 75 




eq. (35) 


P(L,J) L = 1 to 51 


PtiOj) 


Legendre polynomial, see eq. (134) 


J = 1 to 75 






PP(L,J) 


p i%) 


Associated Legendre polynomial, see 


L = 1 to 50 




eq. (135) 


J = 1 to 75 






PMA, PMB 


Pm A /p~N and p mB /p~N 


These are the quantities specified by 
the input as they are more convenient 
than RMA and RMB. 
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FORTRAN Symbol 


Math. Symbol 


Description 


RO 


^ON 


Nuclear radius constant (fermis), see 
eq. (63) 


RC 


#OC 


Charge radius constant (fermis) see 
eq. (66) 


RG 


^OG 


Gaussian radius constant (fermis) see 
eq. (83) 


RHOBC 


PC 


Value of p at which the uniform charge 
density ends, see eq. (74) 


RHOBN 


PN 


Value of p at which the standard po- 
tential falls to half of its initial value, 
see eq. (73) 


RHOBNG 


PG 


Value of p at which the Gaussian ab- 
sorption is centered 


RHOIN(I) 




Input values of p for which the in- 


I = 1 to 250 




tegration interval must change from 
DRHOIN(I-l) to DRHOIN(I). See de- 
scription of subroutine RHOTB) 


ROMAX 




Final value of p in the numerical inte- 
gration 


RHO(I) 


Pi 


Value of p at the i-th interval of inte- 


I = 1 to 250 




gration, see eq. (14) 


RMA, RMB 


pTa A i Pm B 


Values of p at which special form fac- 
tors are matched to standard form fac- 
tors, see eqs. (86) and (87) 


SGMAC(I) 


*c(0i) 


See eq. (136) (square fermis/sterad) 


I = 1 to 75 






SGMAEX(I) 


a ex (^) 


Experimental values of the differential 


I = 1 to 75 




elastic scattering cross section (square 
fermis/sterad) 


SGMATH(I) 


A^l) 


Calculated values of the differential 


I = 1 to 75 




elastic scattering cross section (square 
fermis/sterad), see eq. (34) 


SGMRTH 


°R 


Calculated value of the reaction cross 
section (square fermis) see eq. (132) 


SIGMAO 


^0 


See eqs. (117) and (118) 


SIGMA 1 


01 


See eq. (117) 


SRATIO(I) 


a{0i)/<rc{0i) 


Ratio of calculated to Rutherford cross 


I = 1 to 75 




section 


TA, TV, TW, 




Storage for initial values input for the 


TVS, TWS, TBG, 




parameters 
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FORTRAN Symbol 


Math. Symbol 


Description 


THETAD(I) 


&i 


Scattering angle in center-of-mass sys- 


I = 1 to 75 




tem (degrees) 


THETA(I) 


h 


As above (radians) 


I = 1 to 75 






UCRB(I), UCRM(I) 


U C R(Pi) U CR (p i + %*) 


L-independent part of the real central 


I = 1 to 250 




potential at the beginning and in the 
middle of the i-th interval of integra- 
tion, see eq. (98) 


UCIB(I), UCIM(I) 


Ud(Pi) Uci^ + %*) 


As above for the imaginary central po- 


I = 1 to 250 




tential, see eq. (99) 


USRB(I), USRM(I) 


U S R(Pi) U SR ( Pi + %*) 


As above for the real spin-orbit poten- 


I = 1 to 250 




tial, see eq. (100) 


USIB(I), USIM(I) 


Usi(Pi) u S i(Pi + Ap) 


As above for the imaginary spin-orbit 


I = 1 to 250 




potential, see eq. (101) 


V 


V 


Depth of real central potential (MeV) 


W 


w 


Depth of imaginary central potential 
(MeV) 


VS 


v s 


Real part of spin-orbit potential depth 

(MeV) 


WS 


W S 


Imaginary part of spin-orbit potential 
depth (MeV) 


XC1, XCP1 


xj(p), xj(p) 


Real part of the radial (unnormalized) 
wave function and its first derivative for 
the case L + 1/2 


YC1, YCP1 


yt(p), yj(p) 


As above for the imaginary part and 
the case L + 1/2 


XD1, XDP1 


xj(p), xj(p) 


As above for the real part and the case 
L-l/2 


YD1, YDP1 


vj(p), vjip) 


As above for the imaginary part and 
the case L — 1/2 


XI (L), X1P(L) 


%£ (PmaxJ ? x £ (PmaxJ 


Real part of the radial (unnormalized) 


L = 1 to 51 




wave function and its first derivative for 
the case L + l/2 at the end of a numer- 
ical integration 


Y1(L), Y1P(L) 


Vg {pmax), i/g {Pmax) 


As above for the imaginary part and 


L = 1 to 51 




the case L + 1/2 


X2(L), X2P(L) 


x p (PmaxJ? x g (PmaxJ 


As above for the real part and the case 


L = 1 to 51 




L-l/2 


Y2(L), Y2P(L) 


yg (Pmax)? ilg (Pmax) 


As above for the imaginary part and 


L = 1 to 51 




the case L — 1/2 
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FORTRAN Symbol 


Math. Symbol 


Description 


ZZ 


ZZ' 


Product of the atomic numbers of the 
target nucleus and the incident parti- 
cle. 



VI. Symbolic Listing of the Program 



MAM ROUTINE - SCAT 4 

<HVIVCNA,AR,AI, 

1BR,BI,BG, 

2CHI2 , CHI2P , CHI2PT , CHI2S , CHI2ST , CHI2T , CR1 , CI1 , CR2 , CI2 , 
3DPOLEX, DSGMEX, DRHO, DRHOTN , DRHOL, DV,DW, DA , DVS , DWS, DBG, 
4ECM, ELAB , EPS 1 , EPS2 , EPS3 , EPS4 , ETA , ETA2 , EXSGMR EXSGMI , 
5F , FBAR, FCR, FCI , FFCR, FFCI ,FFCRM, FFCIM , FFSR , FFSI ,FFSRM, FFSIM , 
6FKAY, FMB, FMI , EMU, FN1A , FN2A , FN1B , FN2B , FP , FKAYA, FKAYB, 
7G,GP, 
8HA,HB, 

9rDATA, IFIRST , UN , ILAST , ISPILL 
GOMVONr JMAX, JMAXT, JSPILL , 
1KTRL,KTRLT, 
2L,LMAX,LMAXM, 

3NMAX, NMAXP, NMAXT, NTNPUT , NUMRUN, NUMPRG, NVMAX,NWMAX, NAMAX, NVSMAX, 
4NWSMAX, NV , NW, N A , NVS , NWS , NBGMAX, NBG , 
5P , PP,POLEX,POLTH,PMA,PMB, 

6RC , RO , RHO , RHOBC , RHOBN, RHOIN , RHOMAX, RMA, RMB, RG , RHOBNG, 
7SGMAC, SGMAEX, SGMATH, SGMRTH, SIGMAO , SIGMA1 , SRATIO , 
8THETA,THETAD, TV,TW,TA, TVS ,TWS,TBG, 
9UCRB, UCIB ,UCRM, UCLM , USRB , USIB , USRM, USIM 

oavivosr v,vs, 

1W,WS, 

2X1 , X2 , X1P , X2P , XC1 , XCP1 , XD1 , XDP1 , 

3Y1 , Y2 , YIP , Y2P , YC1 , YCP1 , YD1 , YDP1 , 

4ZZ 

DIMENSION AR(75) ,AI(75) , 
1BR(75) ,BI(75) , 

2CHI2(75) ,CHI2P(75) ,CHI2S(75) ,CR1(51) ,CI1(51) ,CR2(51) ,CI2(51) , 
3DPOLEX(75) ,DSGMEX(75) ,DRHO(250) ,DRHOLN(250) , 
4EXSGMR(51) ,EXSGMI(51) , 

5F(52) ,FBAR(91) ,FCR(75) ,FCI(75) ,FFCR(250) ,FFCI(250) ,FFCRM(250) , 
6FFCIM(250) ,FFSR(250) ,FFSI(250) ,FFSRM(250) ,FFSIM(250) ,FP(51) , 
7G(52) ,GP(51) , 
8IIN(51) , 
9KTRL(13) ,KTRLT(13) 

DIMENSION NUMRUN ( 5 ) , 

1P(51 ,75) ,PP(50,75) ,POLEX(75) ,POLTH(75) , 
2RHO(250) ,RHOLN(250) , 

3SGMAC(75) ,SGMAEX(75) ,SGMATH(75) ,SRATIO(75) , 
4THETA(75) ,THETAD(75) , 

5UCRB(250) ,UCIB(250) ,UCRM(250) ,UCM(250) ,USRB(250) ,USIB(250) , 
6USRM(250) ,USLM(250) , 
7X1(51) ,X2(51) ,X1P(51) ,X2P(51) , 
8Y1(51) ,Y2(51) ,Y1P(51) ,Y2P(51) 

CALL SPILL ( JSPILL , ISPILL , . , . ) 

EPS1= 0.00001 

EPS2= 0.00001 

EPS3= 0.00001 

EPS4 = 0.001 
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READ INPUT TAPE 7 ,10 , (NUMRDN( I )) , 1=1 ,5) 
READ FNPUT TAPE 7,10 ,NUMPRG 
10 PORMAT(I5) 
CALL CTRL4 
GO 
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SUBROUTINE CTRL4 
3 NUMRUN(4)=NUMRUN(4) + 1 
NUMRUN(5) = 
CALL INPT4 
CALL P0T1CH 
35 IF(KTRL(5)) 80,81,80 

80 CALL POP1 

81 CALL SIGZRO 
CALL FSUBC 
CALL EXSGML 

DO 20 NV=1,NVMAX 
IF (NV-1) 102,101,102 
101 V=TV 

GOTO 103 



105 



02 V=V4DV 




03 DO 20 NW= 


1,NWMAX 


IF (NW-1) 


105,104,] 


04 W=TW 




GOTO 109 




05 W=WtDW 




09 DO 20 NA= 


1,NAMAX 


IF (NA-1) 


111 ,110,] 


10 A=TA 




GOTO 112 




11 A=A+DA 




12 DO 20 NVS= 


= 1,NVSMAX 



111 



IF (NVS-1) 114,113,114 

113 VS=TVS 
GOTO 115 

114 VS=VS+DVS 

115 DO 20 NWS=1,NWSMAX 

IF (NWS-1) 117,116,117 

116 WS=TWS 
GOTO 118 

117 WS=WS+DWS 

118 DO 20 NBG=1,NBGMAX 
IF(NBG-l) 120,119,120 

119 BG=TBG 
GOTO 121 

120 BG=BGH)BG 

121 IF (SENSE SWITCH 1) 26,27 
26 REWIND 7 

CALL SAVE(8) 

READ MPUT TAPE 7 ,50 , (LGAR, I = 1 ,6) 

IDATA= NUMRUN(4) 

DO 66 NINPUT=1, ID ATA 

READ MPUT TAPE 7 ,50 , (KTRLT( I ) ,1 = 1,13) 

50 FORMAT (15) 

51 FORMAT (E15.9) 

READ MPUT TAPE 7 ,5 1 , (GAR, I =1 ,27) 

READ MPUT TAPE 7,50, (LGAR, 1 = 1,6), NMAXT 

NT=2*NMAXT-1 

READ MPUT TAPE 7,51, (GAR, I = 1 ,NT) 

READ MPUT TAPE 7,51 ,LGAR 
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IF (KTRLT(5)) 71,70,71 
71 READ INPUT TAPE 7,50, JMAXT 

READ MPUT TAPE 7,51, (GAR, I = 1 ,JMAXT) 
70 IF (KTRLT(2)) 61,66,61 
61 IF(KTRLT(3)) 63,66,63 
63 NT=4*JMAXT 

READ MPUT TAPE 7,51, (GAR, I = 1 ,NT) 
66 CONTINUE 
27 NUMRUN(5)= NUMRUN(5) + 1 

CALL RHOTB 

CALL COULFN 

CALL RMXINC 

CALL PGEN4 

CALL IMTCTR 

CALL CSUBL 

CALL AB 

CALL SGSGCP 

CALL SIGMAR 

IF (KTRL(2)) 33,100,33 
33 CALL CHISQ 
100 CALL OUTPT4 
20 CONTINLE 

GOTO 3 
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SUBROUTINE INPT4 

IF DIVIDE CHECK 100,110 

100 WRITE OUTPUT TAPE 6 , 1 1 

101 PORMAT(59H DIVIDE CHECK TRIGGER FOUND ON AT START OF INPT4 SUBROUT 
1INE) 

CALL LEAVE 

STOP 
110 ISPILL=0 

JSPILL=0 

READ MPUT TAPE 7,10 ,KTRL( 1 ) 

IF (KTRL(l)-lOO) 150,151,151 
151 CALL EXIT 

STOP 
150 READ MPUT TAPE 7 , 10 , (KTRL( I ) ,1=2,13) 
10 FORMAT (15) 

READ MPUT TAPE 7 , 12 ,FMI,FMB,ELAB, ZZ ,RC, V,W,RO, A, VS ,WS,RG,BG, 
1DV,DW, DA , DVS , DWS,DBG 

READ MPUT TAPE 7 , 1 2 , HA , PMA, FN1 A , FN2A , HB , PMB, FN1B , FN2B 

READ MPUT TAPE 7 ,10 ,NVMAX,NVMAX,NAMAX,NVSMAX,NWSMAX,NBGMAX 
12 FORMAT (E15.9) 

TV= V 

TW=W 

TA=A 

TVS=VS 

TWS=WS 

TBG=BG 

READ MPUT TAPE 7,10 ,NMAX 

NMAXP=NMAX-1 

READ MPUT TAPE 7,12, (RHOM( I ) , I = 1 ,NMAX) , (DRHOIN( I ) , I = 1 ,NMAXP) 

C02=FMEFMB 

FMU= (FMI*FMB) / C02 

ECTvT=ELAB* (FMB/C02) 

FKAY= .2 195376 *SQRTF(FMU*ECM) 

T=FKAY*(FMB**. 333333333) 

RHOBN= T*RO 

RHOBNG=T*RG 

RMA=PMA*RHOBN 

RMB=PMB*RHOBN 

RHOBC= T*RC 

ETA= .15805086*ZZ*SQRTF(FMI/ELAB) 

IF DIVIDE CHECK 200,47 

200 WRITE OUTPUT TAPE 6,201 

201 FORMAT(43H MPUT DIVISOR WAS ZERO IN INPT4 SUBROUTINE) 
CALL LEAVE 

STOP 

47 READ MPUT TAPE 7,10 ,LMAXM 
LMAX=LMAXMfl 

DO 147 J = 1,LMAX 
147 IIN(J) = 1 

IF (KTRL(5)) 48,50,48 

48 READ MPUT TAPE 7,10 ,JMAX 

READ MPUT TAPE 7,12, (THETAD( I ) , I = 1 ,JMAX) 
DO 49 I = 1,JMAX 

49 THETA(I)= 0.01745329252*THETAD( I ) 
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50 IF(KTRL(2)) 51,207,51 

51 IF(KTRL(3)) 53,207,53 

53 READ INPUT TAPE 7,12, (SGMAEX( I ) , I=1,JMAX) , (DSGMEX( I ) , I = 1 ,JMAX) 

l(POLEX(I) ,I = 1,JMAX) ,(DPOLEX(I) ,I = 1,JMAX) 
207 IF(ISPILL)202 ,204,202 
202 WRITE OUTPUT TAPE 6,2 03, ISPILL 
203 FORMAT(23H UNDERFLOW OCCURRED AT 15 ,20H IN INPT4 SUBROUTINE) 

204 IF(JSPILL)205 ,210,205 

205 WRITE OUTPUT TAPE 6,2 06, JSPILL 

206 FORMAT(22H OVERFLOW OCCURRED AT 15 ,20H IN INPT4 SUBROUTINE) 
CALL LEAVE 
STOP 
210 RETURN 
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SUBROUTINE P0T1CH 

IF DIVIDE CHECK 30,31 

30 WRITE OUTPUT TAPE 6 ,130 

130 FORMAT (6 OH DIVIDE CHECK TRIGGER FOUND ON AT START OF POT1CH SUBRO 

1UTME) 
CALL LEAVE 
STOP 

31 ISPILL=0 
JSPILL=0 
IKTRETCTRL(13) 
NMAX=NMAX 
NMAXP= NMAX-1 
AMAX=^AMAX-1 

TTA=MAX1F( A , ( (AMAX*DA)+A) ) 
VMAX=NrVMAX-l 

TTV=^V1AX1F(V, ( (VMAX*DV)+V) ) 
WMAXNWMAX-1 

TTW=^IAX1F(W, ( (WMAX*DW)+W) ) 
VSWA5HMVSMAX-1 

TTVS=^LAX1F( VS , ( ( VSMAX*DVS)+VS ) ) 

WSMAXNIWSMAX-1 

TTWS^VIAX1F(WS, ( (WSMAX*DWS)+WS) ) 

BGMAXNBGMAX-1 

TTBG^V1AX1F(BG, ( (BGMAX*DBG)+BG) ) 

FKAYA=TT<AY*TTA 

FKAYB=TKAY*TTBG 

T2^3QRTF ( TTV* * 2 +TIW* * 2 ) /ECM 

T7=TTV/ECM 

T8=TTW/ECM 

IF DIVIDE CHECK 60,61 

60 WRITE OUTPUT TAPE 6,160 

160 FORMAT(26H ECM IS ZERO IN POT1CH SUB) 
CALL LEAVE 
STOP 

61 GOTO (3,3,111 ,15) ,IKTRL 
3 IF(KTRL(l)-2) 24,25,24 

25 IF (RHOIN (NMAX) -RHOBN) 10,10,8 

24 Tl = l./(l.+EXPF( (RHOIN (NMAX) -RHOBN) /FKAYA) ) 

IF DIVIDE CHECK 50,28 
50 WRITE OUTPUT TAPE 6,150 
150 FORMAT(28H FKAYA IS ZERO IN POT1CH SUB) 

CALL LEAVE 

STOP 
28 IF(KTRL(1)-1) 40,41,40 

40 T3= T2*T1 
GOTO 43 

41 T3=T7*T1 

43 IF(T3-EPS4) 42,42,10 

10 WRITE OUTPUT TAPE 6,100, RHOIN (NMAX) , DRHOIN (NMAXP) 

100 FORMAT(13H RHOIN (NMAX) =E 16 . 9 ,2H+ E16.9,46H RHOIN(NMAX) IS TOO SMAL 
1L IN NUCLEAR POTENTIAL) 

RHOIN (NMAX) = RHOIN (NMAX) +DRHOIN (NMAXP) 

GOTO 3 

42 IF(KTRL(1)-1) 8,6,8 
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6 Tl 1= EXPF( - ( (RHOLN (NMAX) -RHOBNG) /FKAYB) * * 2 ) 
IF((T8*T11)-EPS4) 8,8,7 

7 WRITE OUTPUT TAPE 6,103, RHOLN (NMAX) , DRHOLN (MMAXP) 

103 FORMAT(13H RHOLN (NMAX) =E16. 9 ,2H+ E16.9,46H RHOIN (NMAX) IS TOO SMAL 
1L IN NUCLEAR POTENTIAL) 
RHOIN (NMAX) = RHOLN (NMAX) +DRHOLN (NMAXP) 
GOTO 6 

8 GOTO(lll ,15) ,KTRL 
1 1 1 FLMAXLMAXM 

IF(KTRL(l)-2) 29,300,29 
300 IF(FLMAX-(RHOBN+3.)) 12,12,15 
29 T4=l./(l.+EXPF((FLMAX-RHOBN)/FKAYA)) 

IF(KTRL(1)-1) 33,32,33 
33 T5= T2*T4 

GOTO 310 
32 T5=T7*T4 
310 IF(T5-EPS4)13,13,12 

12 WRITE OUTPUT TAPE 6,101 LMAXM 

101 FORMAT (7H LMAXM=I5 ,3H +1,45H LMAXM TOO SMALL BECAUSE OF CENTRAL P 
10TENTIAL) 

LMAX= LMAX+1 
IMAXM= LMAXM+1 
IIN(LMAX) = 1 
GOTO 111 
13 IF(KTRL(1)-1) 17,19,17 

19 T4=EXPF(-((FLMAX-L{HOBNG)/FKAYB)**2) 
IF((T8*T4)-EPS4) 17,17,20 

20 WRITE OUTPUT TAPE 6,200 LMAXM 

200 FORMAT (7H IMAXM=I5 ,3H +1,45H LMAXM TOO SMALL BECAUSE OF CENTRAL P 
10TENTIAL) 
LMAXLMAX+1 
LMAXMJMAXM+1 
IIN(LMAX) = 1 
GOTO 19 

17 T2=SQRTF(TTVS**2+TTWS**2)/ECM 

18 FLMAXLMAXM 

T4= 1 . / ( 1 . +EXPF ( (FLMAX-RHOBN) /FKAYA) ) 
38 T6 = 2.*T2*T4*(FKAYW**2) 
IF(T6-EPS4) 15,15,14 
14 WRITE OUTPUT TAPE 6,102, LMAXM 

102 FORMAT (7H LMAXM=I5 ,3H +1,48H LMAXM TOO SMALL BECAUSE OF SPIN ORB 
1IT POTENTIAL) 

LMAX= LMAX+1 

IMAXM= LMAXM+1 

IIN(LMAX) = 1 

GOTO 18 
15 IF(ISPILL)202 ,204,202 
202 WRITE OUTPUT TAPE 6,2 03, ISPILL 
203 FORMAT(23H UNDERFLOW OCCURRED AT 15 ,14H IN POT1CH SUB) 

204 IF(JSPILL)205 ,210,205 

205 WRITE OUTPUT TAPE 6,2 06, JSPILL 

206 FORMAT(22H OVERFLOW OCCURRED AT 15 ,14H IN POT1CH SUB) 
CALL LEAVE 
STOP 
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210 RETURN 
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SUBROUTINE P0P1 

IF DIVIDE CHECK 1 , 2 

1 WRITE OUTPUT TAPE 6,101 

101 FORMAT (58H DIVIDE CHECK TRIGGER FOUND ON AT START OF POP1 SUBROUT 
1INE) 

CALL LEAVE 
STOP 

2 ISPILL=0 
JSPILL=0 
LMAXPTWAX+1 
DO 20 J = 1,JMAX 

SI2 = l./SINF (THETA( J ) ) 
IF DIVIDE CHECK 3,4 

3 WRITE OUTPUT TAPE 6,103, J 

103 FORMAT (71H DIVISOR SINF THETA IS ZERO IN FIRST DIVISION OF POP1 S 
1UBROUTTNE FOR J=I3) 
CALL LEAVE 
STOP 

4 CO=COSF(THETA(J)) 
P(l ,J) = 1.0 
P(2,J)=CO 

PP(1 ,J) = 0.0 
TWOLPl=3. 
FL=1. 

DO 20 L=1,LMAXP 
TL=FL+1. 

P(L+2,J) = (TWOLPl*CO*P(L+l,J)-FL*P(L,J))/TL 
PP(L+l,J)=TL*SI2*(CO*P(L+l,J)-P(L+2,J)) 
TWOLPl=TWOLPl+ 2 . 
20 FL=TL 

IF (ISPILL) 30,31,30 

30 WRITE OUTPUT TAPE 6,130, ISPILL 

130 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,19H IN POP1 SUBROUTINE) 

31 IF (JSPILL) 32,33,32 

32 WRITE OUTPUT TAPE 6,132, JSPILL 

132 FORMAT (22H OVERFLOW OCCURRED AT 16 ,19H IN POP1 SUBROUTINE) 
CALL LEAVE 
STOP 

33 RETURN 
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SUBROUTINE SIGZRO 
IF DIVIDE CHECK 5 , 6 

5 WRITE OUTPUT TAPE 6,105 

105 FORMAT (6 OH DIVIDE CHECK TRIGGER FOUND ON AT START OF SIGZRO SUBRO 

1UTTNE) 
CALL LEAVE 
STOP 

6 ISPILL = 
JSPILL = 

SIGMA0=-(ETA/ (12.* (ETA* * 2 + 16. )))*(1. + (ETA* *2-48.)/(30.*( (ETA**2 + 16 
1 . ) * * 2 ) ) + ( ETA * * 4 - 1 6 . * ( ETA **2) + 1280.)/(((16.+ ETA **2)**4)*105.)) 
SIGMAO=SIGMAO-ETA+(ETA/ 2 . ) *LOGF(ETA* *2+16.) + ((7./2.) *ATANF(ETA/ 4 . ) 
1 ) - (ATANF (ETA) +ATANF (ETA/ 2 . ) + ATANF (ETA / 3 . ) ) 
SIGMA1=SIGMA0+ATANF(ETA) 
15 IF (ISPILL) 30,31,30 

30 WRITE OUTPUT TAPE 6, 130, ISPILL 
130 FORMAT (23H UNDERFLOW OCCURRED AT 16 ,21H IN SIGZRO SUBROUTINE) 

31 IF (JSPILL) 32,11,32 

32 WRITE OUTPUT TAPE 6, 13 2, JSPILL 

132 FORMAT (22H OVERFLOW OCCURRED AT 16 ,21H IN SIGZRO SUBROUTINE) 

CALL LEAVE 

STOP 
11 RETURN 
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SUBROUTINE FSUBC 

IF DIVIDE CHECK 20,21 

20 WRITE OUTPUT TAPE 6,120 

120 FORMAT (53H DIVIDE TRIGGER FOUND ON AT START OF FSUBC SUBROUTINE) 
CALL LEAVE 
STOP 

21 ISPILL=0 
JSPILL=0 

DO 10 J = 1,JMAX 

SN=(SINF (THETA( J ) / 2 . ) ) * * 2 

FLN=ETA * (LOGF ( SN) ) - 2 . * SIGMA0 

FNO=ETA / ( 2 . * FKAY * ( SN ) ) 

IF DIVIDE CHECK 22,23 

22 WRITE OUTPUT TAPE 6 , 1 2 2 , J 

122 FORMAT (23H DIVISOR IS ZERO FOR J=I3 ,2 OH IN FSUBC SUBROUTINE) 
CALL LEAVE 
STOP 

23 FCR(J) = (-FNO*COSF(FLN)) 
10 FCI(J) = (FNO*SINF(FLN)) 

IF (ISPILL) 24,25,24 

24 WRITE OUTPUT TAPE 6,124, ISPILL 

124 FORMAT (23H UNDERFLOW OCCURRED AT 16 ,20H IN FSUBC SUBROUTINE) 

25 IF (JSPILL) 26,27,26 

26 WRITE OUTPUT TAPE 6,126, JSPILL 

126 FORMAT (22H OVERFLOW OCCURRED AT 16 ,20H IN FSUBC SUBROUTINE) 
CALL LEAVE 
STOP 

27 RETURN 
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SUBROUTINE EXSGML 

IF DIVIDE CHECK 10,11 

10 WRITE OUTPUT TAPE 6 , 1 1 

110 FORMAT (6 OH DIVIDE CHECK TRIGGER POUND ON AT START OF EXSGML SUBRO 
1UTPNE) 
CALL LEAVE 
STOP 

11 ISPILL=0 
JSPILL=0 

1 FL=0. 

EXSGMR( 1 ) =COSF (2.0* SIGMAO ) 

EXSGMI( 1 ) = SINF (2.0* SIGMAO ) 

ETA2=TTTA**2 

ETA2A=2.0*ETA 

LX) 20 L=2,LMAX 

FL=FL+1.0 

TER0=FL**2 

TER1=TER0+ETA2 

TER2= (TER0-ETA2 ) / TER1 

TER3=(ETA2A*FL) /TER.1 

IF DIVIDE CHECK 12,13 

12 WRITE OUTPUT TAPE 6,1 12, L 

112 FORMAT (44H DIVISOR IS ZERO IN EXSGML SUBROUTINE FOR L=I3) 
CALL LEAVE 
STOP 

13 EXSGMR(L) = (TER2*EXSGMR(L-1))-(TER3*EXSGMI(L-1)) 
20 EXSGMI(L) = (TER2*EXSGMI(L-1)) + (TER3*EXSGMR(L-1)) 

IF (ISPILL) 14,15,14 

14 WRITE OUTPUT TAPE 6 , 1 1 4 , ISPILL 

114 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,21H IN EXSGML SUBROUTINE) 

15 IF (JSPILL) 16,17,16 

16 WRITE OUTPUT TAPE 6, 11 6, JSPILL 

116 FORMAT(22H OVERFLOW OCCURRED AT 16 ,21H IN EXSGML SUBROUTINE) 
CALL LEAVE 
STOP 

17 RETURN 
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SUBROUTINE RHOTB 

DRHO(l)=DRHOIN(l) 
RHO(l)=RHOIN(l) 
N=l 
1=1 
20 RHO(I + l)=RHO(I)+DRHOIN(N) 

IF (RHO(I+1)-RHOIN(NMAX))30,50,70 
30 IF(ABSF(RHO(I + l)-RHOIN(N+l))-.5*DRHOIN(N)) 35,35,40 
35 N=XMINOF(N+l,NMAX-l) 

40 DRHO(I + l)=DRHOIN(N) 
1=1+1 
GOTO 20 
50 ILAST=I+1 
60 RHO(ILAST)=RHOIN(NMAX) 

DRHO(ILAST-l)=RHO(ILAST)-RHO(ILAST-l) 
RHOMAX=RHOLN (NMAX) 
DRHOL=DRHOIN (NMAX- 1 ) 
IF(ISPILL) 80,81,80 
80 WRITE OUTPUT TAPE 6,1 80, ISPILL 
180 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,21H IN RHOTB SUBROUTINE) 
81 IF(JSPILL)82,83,82 

82 WRITE OUTPUT TAPE 6,1 82, JSPILL 

182 FORMAT(22H OVERFLOW OCCURRED AT 16 ,21H IN RHOTB SUBROUTINE) 
CALL LEAVE 
STOP 

83 RETURN 

70 IF((RHO(I + 1)-RHOLN(NMAX))-.5*DRHOLN(N))50,50,75 
75 ILAST=I 
GOTO 60 



66 



SUBROUTINE COULFN 

IF DIVIDE CHECK 50,51 

50 WRITE OUTPUT TAPE 6,150 

150 FORMAT (6 OH DIVIDE CHECK TRIGGER FOUND ON AT START OF COULFN SUBRO 
1UTTNE) 
CALL LEAVE 
STOP 

51 ISPILL=0 
JSPILL=0 
IKTRETCTRL(13) 
LMAX=LMAXM+1 
ETA2=ETA**2 
SQ=SQRTF(1.+ETA2) 

1 IJ = 1 
AR(1) = -ETA 
AI(1) = 0. 
AR(2) = -.5*ETA2 
AI(2) = .5*ETA 

2 SI=0. 
SR=0. 

PR= RHOMAX 

DO 10 K=2,49 

T= PR*FLOATF(l-K) 

TR=AR(K) /T 

TI=AI(K)/T 

IF DIVIDE CHECK 52,53 

52 WRITE OUTPUT TAPE 6,152 

152 FORMAT(57H DIVISOR T IS ZERO IN FIRST DIVISION OF COULFN SUBROUTTN 
IE) 





CALL LEAVE 




STOP 


53 


SQN=TR**2+TI**2 




IF(K-2) 4,4,3 


3 


IF(SQN-SQO) 4,4,11 


4 


TR=SR+TR 




TI=SI+TI 




IF(TR-SR) 6,5,6 


5 


IF(TI-SI) 6,13,6 


6 


SR=TR 




SI=TI 




AR(K+1)=0. 




AI(K+1)=0. 




KP=K/2 




DO 7 M=1,KP 




KM=K+1-M 




AR(K+1)=AR(K+1)-AR(M) *AR(KM)+ AI (W) * AI (KM) 




AI(K+1)=AI(K+1)-AI(KM)*AR(M)-AI(M)*AR(KM) 




IF(K-2*KP) 8,9,8 



AR(K+1)=AR(K+1)-.5*(AR(KP+1)**2-AI(KP+1)**2) 

AI(K+1)=AI(K+1)-AR(KP+1)*AI(KP+1) 

FK=.5*FLOATF(K) 

AI(K+1)=AI(K+1)-FK*AR(K) 

AR(K+1)=AR(K+1)+FK*AI(K) 

PR= PR*RHOMAX 
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10 SQO=SQN 
GOTO 101 

11 T=SR**2+SI**2 
IF(T) 105,105,12 

12 IF(ABSF(SQO/T)-EPS3) 13,13,106 

13 GOTO (14 ,15) ,IJ 

1 4 PAR=fflOMA3tETA*LOGF ( 2 . *EHOMAX) 
PHI0R=PAR+SIGMA0-H3R 
PHI0I=SI 

AR(2) = -1.+AR(2) 

13=2 

GOTO 2 

15 PHI1R=PAR+SIGMA1-1.570796325 + SR 
PHI1I=SI 

25 T1=EXPF(-PHI0I) 
T2=EXPF(-PHI1I) 
G(l)=Tl*COSF(PHI0R) 
G(2) = T2*C0SF(PHI1R) 
F1=T1*SINF(PHI0R) 
F2=T2*SINF(PHI1R) 
IF(ABSF(F1*G(2)-F2*G(1)-1./SQ)-EPS1) 31,31,102 

31 IDEC=11 

32 I=LMAX+IDEC 
FBAR( I ) = . 1 
FBAR(I+1) = 0. 
LIMIT=IMAXM+IDEC 
FL=LMAX+11 

T1=SQRTF((FL+1.)**2+ETA2) 
IF(JSPILL) 139,133,139 

139 WRITE OUTPUT TAPE 6, 1390, JSPILL 

1390 FORMAT(23H OVERFLOW2 OCCURRED AT 16 ,21H IN COULFN SUBROUTINE) 
CALL LEAVE 
STOP 
133 DO 33 1 = 1, LIMIT 
L=LMAX+IDEC-I 
FL=L 

T2=SQRTF ( FL* *2 +ETA2 ) 

FBAR(L) = ( (2 . *FL+ 1 .) * (ETA+FL* (FL+ 1 . ) /RHOMAX) *FBAR(L+1)-FL*T1*FBAR(L 
l + 2))/((FL + l.)*T2) 
IF DIVIDE CHECK 54,600 
54 WRITE OUTPUT TAPE 6,154 
154 FORMAT(56H DIVISOR IS ZERO IN SECOND DIVISION OF COULFN SUBROUTINE 

1) 

CALL LEAVE 
STOP 

600 IF(JSPILL) 601,33,601 

601 WRITE OUTPUT TAPE 6,1601 , JSPILL 

1601 FORMAT(22H OVERFLOW OCCURRED AT 16 ,21H IN COULFN SUBROUTINE, 2 4H MU 
1LTIPLY FBAR(I) BY 0.1) 
K=LMAX+IDEC 
FBAR(K)=FBAR(K) * . 1 
JSPILL=0 
GOTO 133 

33 T1=T2 
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ALPHA= 1 . / ( (FBAR( 1 ) *G(2) -FBAR( 2 ) *G( 1 ) ) *SQ) 
IF DIVIDE CHECK 55,43 

55 WRITE OUTPUT TAPE 6,155 

155 FORMAT (55H DIVISOR IS ZERO IN THIRD DIVISION OF COULFN SUBROUTINE 

1) 
CALL LEAVE 

STOP 

43 IMAXP=LMAX+1 
DO 34 I = 1,LMAXP 

34 FBAR( I )=ALPHA*FBAR( I ) 
IF(IDEC-ll) 371,35,371 

371 IF(ABSF(F1/FBAR(1)-1.)-EPS2) 37,37,35 

35 DO 36 I = 1,LMAXP 

36 F(I)=FBAR(I) 
IDEC=IDEC+5 

IF (IDEC-40) 32,32,103 

37 DO 38 I = 1,LMAXP 
IF(ABSF(F(I)/FBAR(I)-1.)-EPS2) 44,44,35 

44 IF DIVIDE CHECK 56,38 

56 WRITE OUTPUT TAPE 6,1 56, L, I 

156 FORMAT(74H DIVISOR FBAR(I)-1. IS ZERO IN FOURTH DIVISION OF COULFN 
1 SUBROUTINE FOR L=I3 ,7H AND 1=13) 

CALL LEAVE 
STOP 

38 CONTINUE 

DO 381 I = 1,MAXP 

381 F(I)=FBAR(I) 

382 T1=SQ 

DO 40 L=1,LMAX 
FL=L 

T2=SQRTF ( ( FL+ 1 . ) * * 2 +ETA2 ) 

G(L + 2) = ((2.*FL+1.)*(ETA+FL*(FL+1.)/RH0WAX)*G(L+1)-(FL+1.)*T1*G(L)) 
1/(FL*T2) 
TS=FL/T1 

IF DIVIDE CHECK 57,45 
57 WRITE OUTPUT TAPE 6,157 

157 FORMAT(58H DIVISOR Tl IS ZERO IN FIFTH DIVISION OF COULFN SUBROUTI 
1NE) 

CALL LEAVE 
STOP 

45 TF(ABSF(F(L)*G(L+1)-F(L + 1)*G(L)-TS)-EPS1) 40,40,104 

40 T1=T2 

41 DO 42 L=1,LMAX 
FL=L 

T=FL**2 

Tl=T/RHOMAX+ETA 

IF DIVIDE CHECK 58,46 
58 WRITE OUTPUT TAPE 6,158 

158 FORMAT (62H DIVISOR RHOMAX IS ZERO IN SIXTH DIVISION OF COULFN SUB 
1ROUTLNE) 

CALL LEAVE 

STOP 

46 T2=SQRTF(T+ETA2) 

FP(L) = (T1*F(L)-T2*F(L + 1))/FL 
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42 GP(L) = (T1*G(L)-T2*G(L + 1))/FL 
IF DIVIDE CHECK 59,47 

59 WRITE OUTPUT TAPE 6,159 

159 FORMAT(60H DIVISOR FL IS ZERO IN SEVENTH DIVISION OF COULFN SUBROU 
1TINE) 

CALL LEAVE 

STOP 
47 IF(ISPILL) 60,61,60 

60 WRITE OUTPUT TAPE 6,1 60, ISPILL 

160 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,21H IN COULFN SUBROUTINE) 

61 IF(JSPILL) 62,63,62 

62 WRITE OUTPUT TAPE 6,1 62, JSPILL 

162 FORMAT(22H OVERFLOW OCCURRED AT 16 ,21H IN COULFN SUBROUTINE) 
CALL LEAVE 
STOP 

63 RETURN 

1 1 WRITE OUTPUT TAPE 6,121 ,RHOMAX,DRHOL 
GOTO (110,110,109,109) ,EKTRL 

109 WRITE OUTPUT TAPE 6,114 
GOTO 13 

102 WRITE OUTPUT TAPE 6,122 ,RHOMAX,DRHOL 
GO TO( 110,110,111,111) ,IKTRL 

111 WRITE OUTPUT TAPE 6,114 
GOTO 31 

103 WRITE OUTPUT TAPE 6 , 123 ,RHOMAX,DRHOL 
GOTO (110,110,112 ,112) ,EKTRL 

112 WRITE OUTPUT TAPE 6,114 
GOTO 382 

104 WRITE OUTPUT TAPE 6,124 ,RHOMAX,DRHOL ,L 
GOTO (110,110,113,113) ,EKTRL 

113 WRITE OUTPUT TAPE 6,114 
GOTO 40 

105 WRITE OUTPUT TAPE 6,125 ,RHOMAX,DRHOL 
GOTO (110,110,115 ,115) ,EKTRL 

115 WRITE OUTPUT TAPE 6,114 
GOTO 12 

106 WRITE OUTPUT TAPE 6,126 ,RHOMAX,DRHOL 
GOTO (110,110,116,116) ,EKTRL 

116 WRITE OUTPUT TAPE 6,114 
GOTO 13 

110 RHOMAX=FiHOMAXtDRHOL 
GOTO 1 
121 FORMAT(18H INCREASE RHO MAX=E11.4 ,2H+ Ell. 4 ,35H A OR B SERIES CONV 
1ERGES TOO SLOWLY) 
122 E0RMAT(18H ENCREASE RHO MAX=E11.4 ,2H+ E11.4 ,22H BAD INITIAL WRONSK 
HAN) 

123 FORMAT(18H LNCREASE RHO MAX=E11.4 ,2H+ Ell. 4 ,24H L TOO LARGE IN FBA 
1R (L)) 

124 FORMAT(18H LNCREASE RHO MAX=E11.4 ,2H+ Ell. 4 ,21H BAD WRONSKIAN FOR 
1L=I3) 

125 PORMAT(67H SERIES IN PHI0 OR PHI1 IS ZERO, CHECK DATA, IF OK LNCRE 
1ASE RHOMAX=El 1 . 4 , 2H+ El 1 . 4 ) 

126 PORMAT(52H A OR B SERIES DIVERGES TOO QUICKLY LNCREASE RHOMAX=Ell . 
14,2H+ Ell. 4) 
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1 1 4 FORMAT( 4 2H RHOMAX INCREASE NOT PERMITTED BY KTRL (13)) 



SUBROUTINE RMXLNC 

3 IF (RHOMAX-RHO( ILAST ) ) 1,2,1 

1 ILAST=ILAST+1 
RHO(ILAST)=RHO(ILAST-l)+DRHOL 
DRHO ( ILAST - 1 ) =DRHOL 

GOTO 3 

2 RETURN 
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SUBROUTINE PGEN4 

IF DIVIDE CHECK 60,61 

60 WRITE OUTPUT TAPE 6,160 

160 FORMAT (59H DIVIDE CHECK TRIGGER POUND ON AT START OF PGEN4 SUBROU 
1TINE) 
CALL LEAVE 
STOP 

61 ISPILL=0 
JSPILL=0 
IF(KTRL(1)) 3,4,3 

3 KTRL(7) = 
KTRL(8) = 
KTRL(9) = 
KTRL(10) = 

4 T1=V/ECM 
T2=\V/ECM 
T10=VS/ECM 
T11=WS/ECM 
T12=FKAY*BG 
T3 = 2.*FKAY/A 

IF DIVIDE CHECK 62,65 

62 WRITE OUTPUT TAPE 6,162 

162 FORMAT (65H DIVISORS ECM OR A WERE WRONGLY MPUT AS ZERO IN PGEN4 
1SUBROUTINE) 

CALL LEAVE 
STOP 
65 T4=T10*T3 
T5=T11*T3 
T6=FKAY*A 
T7=ETA/RHOBC 
IF DIVIDE CHECK 63,64 

63 WRITE OUTPUT TAPE 6,163 

163 FORMAT(61H DIVISOR RHOBC IS ZERO IN SECOND DIVISION OF PGEN4 SUBRO 
1UTTNE) 

CALL LEAVE 
STOP 

64 T8=RHOBC**2 
T9=ETA*2. 
1=1 

40 EX=EXPF( (RHO( I)-RHOBN) /T6) 

IF DIVIDE CHECK 80,66 
80 WRITE OUTPUT TAPE 6,165 
165 FORMAT (58H QUANTITY T6 IS ZERO IN THIRD DIVISION OF PGEN4 SUBROUT 
1INE) 





CALL LEAVE 




STOP 


66 


K=l 


41 


IF(I-l) 42,43,42 


42 


IF(DRHO(I)-DRHO(I-l)) 43, 


43 


HDRHO=DRHO( I ) * . 5 




DEX=EXPF (HDRHO/T6 ) 


44 


IF(KTRL(l)-2)53,52,53 


52 


IF(RHO(I)-RHOBN) 54,55,55 


54 


Sl = 1.0 
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GOTO 68 

55 S1=0.0 

GOTO 68 

53 S1 = 1./(1.+EX) 

IF DIVIDE CHECK 67,68 

67 WRITE OUTPUT TAPE 6,167 

167 FORMAT(60H DIVISOR l.+EX IS ZERO IN FOURTH DIVISION OF PGEN4 SUBRO 
1UTME) 
CALL LEAVE 
STOP 

68 S2=EX*(S1**2) 
S4=S2/RHO(I) 

IF DIVIDE CHECK 69,70 

69 WRITE OUTPUT TAPE 6,169,1 

169 FORMAT(58H DIVISOR RHO IS ZERO IN FIFTH DIVISION OF PGEN4 SUBROUTI 
1NE) 

CALL LEAVE 
STOP 

70 IF (RHO(I)-RHOBC) 9,9,10 

9 S3=T7*(3.-(RHO(I)**2)/T8) 
GOTO 11 

10 S3=T9/RHO(I) 

11 IF (KTRL(7)) 350,300,350 

300 UCRB(I)=-1.-T1*S1+S3 
FFCR(I) = S1 

301 IF (KTRL(8)) 355,302,355 
302 IF(KTRL(1)-1) 309,308,309 

308 Sl=EXPF(-((RHO(I)-RHOBNG)/T12)**2) 
IF DIVIDE CHECK 82,309 

82 WRITE OUTPUT TAPE 6,182 

182 FORMAT(22H BG IS ZERO IN PGEN SR) 

CALL LEAVE 

STOP 

309 UCIB(I)=-T2*S1 
FFCI(I) = S1 

303 IF (KTRL(9)) 360,304,360 

304 USRB(I)=T4*S4 
FFSR(I) = S4 

305 IF (KTRL(ll)) 501,500,501 
500 IF (KTRL(10))365 ,306,365 

306 USIB(I)=T5*S4 
FFSI(I) = S4 

307 IF (I-ILAST) 50,200,200 

350 ITT=1 
GOTO 340 

355 ITT=2 

GOTO 340 
340 ITQ=1 

IF(ITT-l) 380,380,381 
380 IF (KTRL(7)-1) 352,351,352 

351 TW=^IA 
TRVLRMA 
TN1=FN1A 
TN2=FN2A 
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GOTO 400 


352 


TH=HB 




TRM^MB 




TN1=FN1B 




TN2=FN2B 




GOTO 400 


381 


IF (KTRL(8)-1) 352,351,352 


400 


IF (RHO(I)-RHOBN) 410,410,411 


410 


TTR=TN1 




GOTO 412 


411 


TTN=TN2 


412 


T20^HO(I)/RHOBN 




IF (TTN*LOGF(T20)-80.) 403,403,409 


403 


TQ= (T20 * *TTN- 1 . ) *RHOBN/ (TTN*FKAY*A) 




IF DIVIDE CHECK 405,406 


405 


TG=T20 * * (RHOBN/ (FKAY*A) ) 




GOTO 407 


406 


IF (TQ-80.) 408,408,409 


408 


TG=EXPF(TQ) 




GOTO 407 


409 


TF=0. 




GOTO 422 


407 


TFN=1./(1.+TG) 




IF (RHO(I)-TRM) 420,420,419 


419 


TF=TFN 




GOTO 418 


420 


T21^iHO(I)/TRM 




THH=TH*(1. + (2.*T21))*((1. -T21)**2) 




TF=TFN* ( 1 . +THH) 


418 


TFF=TF 


421 


GOTO (422 ,423) ,ITQ 


422 


GOTO (425 ,426,427,428) , ITT 


425 


FFCR(I)=TF 




UCRB( I) = -1.-T1*FFCR( I) + S3 



GOTO 301 
426 FFCI(I)=TF 

UCIB(I)=-T2*FFCI(I) 

GOTO 303 
427 FFSR(I)=TF 

IF (ITQ-1) 470,470,471 
471 USRB(I)=FKAY*A*T4*FFSR(I) 

GOTO 305 
470 USRB(I)=(T4/2.)*FFSR(I) 

GOTO 305 
428 FFSI(I)=TF 

IF (ITQ-1) 472,472,473 
473 USIB(I)=FKAY*A*T5*FFSI(I) 

GOTO 307 
360 ITT=3 

IF (KTRL(9)-1) 431,431,430 
430 ITQ=1 

GOTO 352 
365 ITT=4 

IF (KTRL(lO)-l) 431,431,430 
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472 USIB(I) = (T5/2.)*FFSI(I) 

GOTO 307 
431 ITQ=2 

GOTO 351 
423 T23=(RHOBN/ (FKAY*A) ) * ( T20**TTN) *TG* ( (TFN/RHO( I ) ) * * 2 ) 

T25=T23 

IF(RHO(I)-TRM) 460,460,461 

460 T24=6.*TH*(1.-T21)/(TRM**2) 
T25 = (T24*TFN) + ((1.+THH)*T23) 

461 TF=T25 

IF(ITT-3) 427,427,428 
501 T30 = 0.004927*ETA*ECM 

IF(RHO(I)-RHOBC) 502,502,503 

502 SOCOUL=T30/(RHOBC**3) 
GOTO 504 

503 SOCOUL=T30/(RHO(I)**3) 

504 USRB(I)=USRB(I)+SOCOUL 
GOTO 500 

50 1=1+1 

EX=EX*DEX 

RHOVL4TEIO( I -l)+HDRHO 
IF(KTRL(l)-2) 153,152,153 

152 IF(RHOVHIHOBN)34,35,35 

34 Sl = 1.0 
GOTO 72 

35 S1=0.0 
GOTO 72 

153 S1 = 1./(1.+EX) 

IF DIVIDE CHECK 71,72 

71 WRUE OUTPUT TAPE 6,171 

171 FORMAT(54H DIVISOR 15 ZERO IN SIXTH DIVISION OF PGEN4 SUBROUTINE) 
CALL LEAVE 
STOP 

72 S2=EX*(S1**2) 
S4=S2/RHOM 

IF DIVIDE CHECK 73,74 

73 WRITE OUTPUT TAPE 6,173 

173 FORMAT (62H QUANTITY RHOM IS ZERO IN SEVENTH DIVISION OF PGEN4 SUB 
1ROUTPNE) 
CALL LEAVE 
STOP 

74 IF (RHO^LRHOBC) 21,21,22 

21 S3=T7*(3.-(RHOM**2)/T8) 
GOTO 23 

22 S3=T9/RHOM 

23 IF (KTRL(7))1350, 1300, 1350 

1300 UCRM(I-1)=-1.-T1*S1+S3 
FFCRM(I-1)=51 

1301 IF (KTRL(8)) 1355,1302,1355 

1302 IF(KTRL(1)-1) 1309,1308,1309 

1308 S 1=EXPF( - ( (RHCVLRHOBNG) / T12 ) * * 2 ) 

1309 UOM(I-l)=-T2*Sl 
FFCIM(I-1) = S1 

1303 IF (KTRL(9)) 1360,1304,1360 
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1304 USRM(I-1)=T4*S4 
FFSRM(I-1) = S4 

1305 IF (KTRL(Il)) 1501,1500,1501 
1500 IF (KTRL(I0))1365 ,1306,1365 

1306 USIM(I-1)=T5*S4 
FFSIM(I-1) = S4 

1307 IF (K-10) 24,40,40 

1350 ITT=1 
GOTO 1340 

1355 ITT=2 

GOTO 1340 
1340 ITQ=1 

IF (ITT-1)1380, 1380, 1381 

1380 IF (KTRL(7)-1) 1352,1351,1352 

1351 TH=^IA 
THVHMA 
TN1=FN1A 
TN2=FN2A 
GOTO 1400 

1352 TH=HB 
TKVHiMB 
TN1=FN1B 
TN2=FN2B 
GOTO 1400 

1381 IF (KTRL(8)-1) 1352,1351,1352 
1400 IF (RHOM-RHOBN) 1410,1410,1411 

1410 TTT^=TN1 
GOTO 1412 

1411 TTO=TN2 

1412 T20^HOM/RHOBN 

IF (TTN*LOGF(T20)-80.) 1403,1403,1409 
1403 TQ=(T20**TTN- 1 . ) *RHOBN/ (TTN*FKAY*A) 
IF DIVIDE CHECK 1405,1406 

1405 TG=T20**(RHOBN/(FKAY*A)) 
GOTO 1407 

1406 IF (TQ-80.) 1408,1408,1409 

1408 TG=EXPF(TQ) 
GOTO 1407 

1409 TF=0. 
GOTO 1422 

1407 TFN=1./(1.+TG) 

IF (RHOM-TRM) 1420,1420,1419 

1419 TF=TFN 
GOTO 1418 

1420 T21^iHOM/TRM 

TRH=TH*(1. + (2.*T21))*((1. -T21)**2) 
TF=TFN* ( 1 . +THH) 
1418 TFF=TF 

1421 GOTO (1422 ,1423) ,ITQ 

1422 GOTO (1425 ,1426,1427,1428) , ITT 

1425 FFCRM(I-1)=TF 
UCRM(I-1)=-1.-T1*FFCRM(I-1)+S3 
GOTO 1301 

1426 FFCIM(I-1)=TF 
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UCIM(I-1)=-T2*FFCIM(I-1) 
GOTO 1303 

1427 FF3RM(I-1)=TF 

IF (ITQ-1) 1470,1470,1471 

1471 USRM(I-1)=FKAY*A*T4*FF3RM(I-1) 
GOTO 1305 

1470 USRM(I-1) = (T4/2.)*FFSRM(I-1) 
GOTO 1305 

1428 FFSIM(I-1)=TF 

IF (ITQ-1) 1472,1472,1473 
1473 USIM(I-1)=FKAY*A*T5*FFSIM(I-1) 

GOTO 1307 
1360 ITT=3 

IF (KTRL(9)-1) 1431,1431,1430 

1430 ITQ=1 
GOTO 1352 

1365 IIT=4 

IF (KTRL(lO)-l) 1431,1431,1430 

1472 USIM(I-1) = (T5/2.)*FFSIM(I-1) 
GOTO 1307 

1431 ITQ=2 
GOTO 1351 

1423 T23=(RHOBN/(FKAY*A) ) * ( T20**TTN) *TG* ( (TFN/RHOM) **2) 
T25=T23 
IF(RHOM-TRM) 1460,1460,1461 

1460 T24=6.*TH*(1.-T21)/(TRM**2) 
T25=(T24*TFN) + ( ( 1 . +THH) * T23 ) 

1461 TF=T25 

IF (ITT-3) 1427,1427,1428 

1501 T30 = 0.004927*ETA*ECM 

IF (RHOM-RHOBC) 1502,1502,1503 

1502 SOCOUI^T30/(RHOBC**3) 
GOTO 1504 

1503 SOCOUL-T30/(RHOM**3) 

1504 USRM(I-l)=USRM(I-l)+SOCOUL 
GOTO 1500 

24 K=K+1 

EX=EX*DEX 
GOTO 42 
200 IF(ISPILL) 75,76,75 

75 WRITE OUTPUT TAPE 6,175,ISPILL 

175 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,2 OH IN PGEN4 SUBROUTINE) 

76 IF (JSPILL) 77,51,77 

77 WRITE OUTPUT TAPE 6,177, JSPILL 

177 FORMAT(22H OVERFLOW OCCURRED AT 16 ,2 OH IN PGEN4 SUBROUTINE) 
CALL LEAVE 
STOP 
51 RETURN 
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SUBROUTINE INTCTR 

D01 L=1,LMAX 

IFIRST=IIN(L) 

T=^HO(IFIRST)**(L-l) 

XCl=T*RHO( IFIRST ) 

XD1=XC1 

FI^=L 

XCP1=FL*T 

XDP1=XCP1 

YC1=0. 

YD1=0. 

YCP1 = 0. 

YDP1 = 0. 

CALL RKINT 

X1(L)=XC1 

X2(L)=XD1 

Y1(L)=YC1 

Y2(L)=YD1 

X1P(L)=XCP1 

X2P(L)=XDP1 

Y1P(L)=YCP1 

Y2P(L)=YDP1 

REIURN 
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SUBROUTINE RMNT 

IF DIVIDE CHECK 10,11 
10 WRITE OUTPUT TAPE 6,1 10, L, I 
110 PORMAT(66H DIVIDE CHECK TRIGGER FOUND ON AT START OF RKTNT SUBROUT 
1INE FOR L=I3 ,7H AND 1=13) 
CALL LEAVE 
STOP 

11 ISPILL=0 
JSPILL=0 

1 FL=L-1 

F2L=-1.-FL 
F3L=FL*(FL+1.) 

TB=UCRB(IFIRST)+F3L/(RHO(IFIRST)**2) 
IF DIVIDE CHECK 12,13 

12 WRITE OUTPUT TAPE 6,1 12, L, I 

112 FORMAT(76H DIVISOR RHO(IFIRST) **2 IS ZERO IN FIRST DIVISION OF RKI 
INT SUBROUTINE FOR L=I3 ,7H AND 1=13) 
CALL LEAVE 
STOP 

13 PCB=TB+USRB(IFIRST)*FL 
PDB=TB+USRB ( IFIRST ) * F2L 
QCB=UCIB(IFIRST)+USIB(IFIRST)*FL 
QDB=UCIB ( IFIRST)+USIB ( IFIRST ) * F2L 
IK=ILAST-1 

DO 6 I=IFIRST,IK 
2 BDRHO=.5*DRHO(I) 

DRH02=(DRHO( I ) * * 2) * . 5 
RHOM=RHO( I ) +HDRHO 
TM=UCRM( I ) +F3L / (RHOM* * 2 ) 
IF DIVIDE CHECK 14,15 

14 WRITE OUTPUT TAPE 6,1 14, L, I 

114 FORMAT(70H DIVISOR RHOM** 2 IS ZERO IN SECOND DIVISION OF RKINT SUB 
1ROUTTNE FOR L=I3 , 7H AND 1=13) 
CALL LEAVE 
STOP 

15 PCM=TM+USRM(I)*FL 
PDM=TM+USRM( I ) * F2L 
QCM=UCIM ( I ) +USIM ( I ) * FL 
QDM=UCIM ( I ) +USIM ( I ) * F2L 
XCPP1=PCB*XC1-QCB*YC1 
YCPP1=QCB*XC1+PCB*YC1 
XDPP1=PDB*XD1-QDB*YD1 
YDPP1=QDB*XD1 T PDB*YD1 
XC2=XC1+XCP1 *HDRHO 
YC2=YC1+YCP1 *HDRHO 
XD2=XD1+XDP1 *HDRHO 
YD2=YD1+YDP1 *HDRHO 
XCPP2=PCM*XC2-QCM*YC2 
YCPP2=QCM*XC2^PCM*YC2 
XDPP2=PDM*XD2-QDM*YD2 
YDPP2=QDM*XD24PDM*YD2 
DRH04=.5*DRH02 
SDRHO= . 3 3 3 3 3 3 3 3 *HDRHO 
XC3=XC2+XCPP1 *DRH04 
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YC3=YC2+YCPP1 *DRH04 
XD3=XD2+XDPP1 *DRH04 
YD3=YD2+YDPP1 *DRH04 
XCPP3=PCM*XC3-QCW*YC3 
YCPP3=QCM*XC3+PCM*YC3 
XDPP3=PDM*XD3-QDM*YD3 
YDPP3=QDM*XD3-ffDW*YD3 
XC4=XC2+XCPP2 *DRH02+XCP1 *HDRHO 
YC4=YC2+YCPP2*DRH02+YCP1 *HDRHO 
XD4=XD2+XDPP2 *DRH02+XDP1 *HDRHO 
YD4=YD2+YDPP2 *DRH02+YDP1 *HDRHO 
TB=UCRB( I + 1)+F3L/ (RHO( I + 1)**2) 
IF DIVIDE CHECK 16,17 

16 WRITE OUTPUT TAPE 6,1 16, L, I 

116 PORMAT(74H DIVISOR RHO(I + l)**2 IS ZERO IN THIRD DIVISION FOR RKTNT 
1 SUBROUTINE FOR L=I3 ,7H AND 1=13) 
CALL LEAVE 
STOP 

17 PCB=TB+USRB(I + 1)*FL 
PDB=TB+USRB ( I + 1) *F2L 
QCB=UCIB(I+1)+USIB(I + 1)*FL 
QDB=UCI8 ( I +1HUSIB ( I + 1 ) *F2L 
XCPP4=PCB*XC4-QCB*YC4 
YCPP4=QCB*XC4+PCB*YC4 
XDPP4=PDB*XD4-QDB*YD4 
YDPP4=QDB*XD4+PDB*YD4 
SXO=XCPP2+XCPP3 
SYC=YCPP2+YCPP3 
SXD=XDPP2+XDPP3 
SYD=YDPP2+YDPP3 
TXOSXC+XCPP1 
TYCHSYGfYCPPl 
TXD=SXD+XDPP1 
TYD=SYD+YDPP1 

TXCl=XCl-fDRHO( I ) * (XCP1+SDRH0*TXC) 
TYCl=YCl-fDRHO( I ) * (YCP1+SDRH0*TYC) 
TXDl=XDl-fDRHO( I ) * (XDP1+SDRH0*TXD) 
TYD1=YD14DRW0( I ) * (YDP1+SDRH0*TYD) 
TXCP1=XCPHSDRH0* (TXC+-SXC+XCPP4) 
TYCP1=YCP1+SDRH0* (TYG+-SYC+YCPP4) 
TXDP1=XDP1+SDRH0* (TXD+-SXD+XDPP4) 
TYDP1=YDP1+SDRH0* (TYDfSYD+YDPP4) 
IF (JSPILL) 20,21,20 

20 REN0RVHV[AX1F(ABSF(XC1) ,ABSF(YC1) ,ABSF(XCP1) ,ABSF(YCP1) ,ABSF(XD1) , 
lABSF(YDl) ,ABSF(XDP1) ,ABSF(YDP1)) 
XC1=XC1 /RENORM 
YC1=YC1/RFN0RM 
XCP1=XCP1 /RENORM 
YCP1=YCP1 /RENORM 
XD1=XD1/RFN0RM 
YD1=YD1/RFN0RM 
XDP1=XDP1/REN0RM 
YDP1=YDP1/REN0RM 
WRITE OUTPUT TAPE 6 ,200 ,RENORM,L ,RHO( I ) 
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200 FORMAT(24H RENORMALEZATION FACTOR=E16.9 ,22H IN RKINT FOR CODED L=I 
13, 9H AND RHO=E16.9) 

JSPILL=0 

GO T02 
21 XC1=TXC1 

YC1=TYC1 

XD1=TXD1 

YD1=TYD1 

XCP1=TXCP1 

YCP1=TYCP1 

XDP1=TXDP1 

YDP1=TYDP1 
6 CONTINUE 

IF (ISPILL) 30,31,30 

30 WRITE OUTPUT TAPE 6,130, ISPILL, L, I 

130 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,2 7H IN RKINT SUBROUTINE FOR L= 
113 .7H AND 1=13) 

31 IF (JSPILL) 32,4,32 

32 WRITE OUTPUT TAPE 6,132, JSPILL, L, I 

132 FORMAT(22H OVERFLOW OCCURRED AT 16 ,27H IN RKTNT SUBROUTINE FOR L=I 
13, 7H AND 1=13) 
CALL LEAVE 
STOP 
4 RETURN 
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SUBROUTINE CSUBL 

IF DIVIDE CHECK 50,51 

50 WRITE OUTPUT TAPE 6,150 

150 FORMAT (59H DIVIDE CHECK TRIGGER POUND ON AT START OF CSUBL SUBROU 
1TINE) 
CALL LEAVE 
STOP 

51 ISPILL=0 
JSPILL=0 

DO 40 L=1,LMAX 

XN0RM1=MAX1F(ABSF(X1(L))*ABSF(Y1(L)) ,ABSF(X1P(L) ) , ABSF(Y1P(L) ) ) 

TX1L=N1 ( L ) /XNORM1 

TY1L=Y1 ( L ) /XNORM1 

TX1PL=N1P ( L ) /XNORM1 

TY1PL=Y1P ( L ) /XNORM1 

FN0R]VHV1AX1F(F(L) ,G(L) ,FP(L) ,GP(L)) 

TFL=F(L)/FNORM 

TGL=G(L) /FNORM 

TFPL=FP ( L ) /FNORM 

TGPL=GP ( L ) /FNORM 

C01=TFL*TY1PL-TFPL*TY1L 

C02=TFPL*TX1L-TFL*TX1PL 

C03=TYlL*TGPL-TYlPL*TGLfTXlL*TFPL-TXlPL*TFL 

C04=TXlPL*TGL-TXlL*TGPLfTYlL*TFPL-TYlPL*TFL 

CO7=1.0/(CO3**2+CO4**2) 

IF DIVIDE CHECK 52,53 

52 WRITE OUTPUT TAPE 6,152 

152 FORMAT(54H DIVISOR IS ZERO IN FIRST DIVISION OF CSUBL SUBROUTINE) 
CALL LEAVE 
STOP 

53 CR1(L) = (C01*C03+C02*C04)*C07 
CI1(L) = (C02*C03-C01*C04)*C07 

XN0RM2=MAX1F(ABSF(X2(L)) ,ABSF(Y2(L)) ,ABSF(X2P(L) ) , ABSF(Y2P(L) ) ) 
TX2L=N2 ( L ) /XNORM2 

TY2L=Y2 ( L ) /XNORM2 

TX2PL=N2P ( L ) /XNORM2 

TY2PL=Y2P ( L ) /XNORM2 

C01=TFL*TY2PL-TFPL*TY2L 

C02=TFPL*TX2L-TFL*TX2PL 

C03=TY2L*TGPL-TY2PL*TGLfTX2L*TFPL-TX2PL*TFL 

C04=TX2PL*TGL-TX2L*TGPLfTY2L*TFPL-TY2PL*TFL 

C07= 1.0/ (C03**2+C04* * 2 ) 

IF DIVIDE CHECK 54,55 

54 WRITE OUTPUT TAPE 6,154 

154 FORMAT (55H DIVISOR IS ZERO IN SECOND DIVISION OF CSUBL SUBROUTINE 

1) 

CALL LEAVE 
STOP 

55 CR2(L) = (C01*C03+C02*C04)*C07 
40 CI2(L) = (C02*C03-C01*C04)*C07 
IF (ISPILL) 56,57,56 
56 WRITE OUTPUT TAPE 6, 15 6, ISPILL ,L 

156 FORMAT (23H UNDERFLOW OCCURRED AT 16 ,27H IN CSUBL SUBROUTINE FOR L 
1=13) 
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57 IF (JSPILL) 58,59,58 

58 WRITE OUTPUT TAPE 6,158, JSPILL , L 

158 FORMAT (22H OVERFLOW OCCURRED AT 16 ,2 7H IN CSUBL SUBROUTINE FOR L= 
113) 

CALL LEAVE 
STOP 

59 RETURN 
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SUBROUTINE AB 

IF DIVIDE CHECK 1 , 2 

1 WRITE OUTPUT TAPE 6,101 

101 FORMAT (56H DIVIDE CHECK TRIGGER POUND ON AT START OF AB SUBROUTTN 
IE) 

CALL LEAVE 
STOP 

2 ISPILL=0 
JSPILL=0 
FKAYD=1./FKAY 

IF DIVIDE CHECK 3,4 

3 WRITE OUTPUT TAPE 6,103 

103 FORMAT(38H DIVISOR FKAY IS ZERO IN AB SUBROUTINE) 
CALL LEAVE 
STOP 

4 DO 20 J = 1,JMAX 
ASUMR=0. 
ASUMI=0. 
BSUMR=0. 
BSUMI=0. 

DO 10 L=1,LMAX 

FL=L 

ATR1=FL*CR1 (L) + (FL- 1 . ) *CR2 (L) 

ATI1=FL*CI1(L) + (FL-1.)*CI2(L) 

BTR1=CR1 ( L) -CR2 ( L ) 

BTI1=CI1(L)-CI2(L) 

ATR2=ATR1*EXSGMR(L) - (ATI1*EXSGMI(L) ) 

ATI2=ATR1*EXSGMI( L) + ( ATI1 *EXSGMR( L ) ) 

BTR2=BTR1*EXSGMR(L) - (BTI1*EXSGMI(L) ) 

BTI2=BTR1 *EXSGMI ( L ) + ( BTI 1 *EXSGMR( L ) ) 

ASUMR=^SUMR+ (ATR2*P ( L , J ) ) 

ASUMI=ASUMI+ ( ATI2 *P ( L , J ) ) 

BSUMR=P3SUMR+ (BTR2*PP ( L , J ) ) 
1 BSUMI=BSUMI+ (BTI2 *PP ( L , J ) ) 

AR(J)= FCR(J) + (FKAYD*ASUMR) 

AI ( J ) =FCI ( J ) + (FKAYD* ASUMI) 

BR(J)= FKAYD*BSUMI 
20 BI(J)= -FKAYD*BSUMR 

IF (ISPILL) 30,31,30 
30 WRTTE OUTPUT TAPE 6,130, ISPILL 
130 FORMAT(23H UNDERFLOW OCCURRED AT 16 ,17H IN AB SUBROUTINE) 
31 IF (JSPILL) 32,33,32 

32 WRITE OUTPUT TAPE 6, 13 2, JSPILL 

132 FORMAT (22H OVERFLOW OCCURRED AT 16 ,17H IN AB SUBROUTINE) 
CALL LEAVE 
STOP 

33 RETURN 
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SUBROUTINE SGSGCP 

IF DIVIDE CHECK 10,11 

10 WRITE OUTPUT TAPE 6 , 1 1 

110 FORMAT (6 OH DIVIDE CHECK TRIGGER POUND ON AT START OF SGSGCP SUBRO 
1UTTNE) 
CALL LEAVE 
STOP 

11 ISPILL=0 
JSPILL=0 

DO 5 J = 1,JMAX 

SGMA r ffl(J)=AR(J)**2. + AI(J)**2.+BR(J)**2. + BI(J)**2. 
POLTH( J)= ( 2 . * ( AR( J ) *BR( J)+AI ( J ) * BI ( J ) ) ) /SGMATH( J ) 
IF DIVIDE CHECK 12,13 

12 WRITE OUTPUT TAPE 6 , 1 1 2 , J 

112 FORMAT(30H DIVISOR SGMATH IS ZERO FOR J = 13,21H IN SGSGCP SUBROUTTN 
IE) 

CALL LEAVE 
STOP 

13 SGMAC(J)=FCR(J)**2. + FCI(J)**2. 
IF(ETA) 7,7,8 

8 SRATIO(J)=SGMATH(J)/SGMAC(J) 
IF DIVIDE CHECK 14,15 

14 WRITE OUTPUT TAPE 6 , 1 1 4 , J 

114 FORMAT(29H DIVISOR SGMAC IS ZERO FOR J = 13,21H IN SGSGCP SUBROUTINE 

1) 

CALL LEAVE 
STOP 

15 GOTO 5 

7 SRATIO(J) = 0. 
5 CONTINUE 

IF (ISPILL) 16,17,16 

16 WRITE OUTPUT TAPE 6, 11 6, ISPILL 

116 FORMAT (23H UNDERFLOW OCCURRED AT 16,21H IN SGSGCP SUBROUTINE) 

17 IF (JSPILL) 18,19,18 

18 WRITE OUTPUT TAPE 6, 11 8, JSPILL 

118 FORMAT(22H OVERFLOW OCCURRED AT 16,21H IN SGSGCP SUBROUTINE) 
CALL LEAVE 
STOP 

19 RETURN 
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SUBROUTINE SIGMAR 
ISPILL=0 
JSPILL=0 
FL=0. 
SGMRTH=0. 

CPI=( 12. 56637060) /(FKAY**2) 
DO 20 L=I ,LMAX 

SGMKHtSGMRTIPrFL*(C12 (L) -(C12 (L))**2 - (CR2(L) ) **2) 
FL=FL+1.0 
20 SGMRTItSGMRTIFrFL*(CIl(L)-(CIl(L))**2-(CRl(L))**2) 
SGMRTEfcCPI *SGMRTH 
IF(ISPILL) 10,11,10 

10 WRITE OUTPUT TAPE 6,1 10, ISPILL 

110 FORMAT(23H UNDERFLOW OCCURRED AT 16,21H IN SIGMAR SUBROUTINE) 

11 IF(JSPILL) 12,13,12 

12 WRITE OUTPUT TAPE 6,1 12, JSPILL 

112 FORMAT(22H OVERFLOW OCCURRED AT 16,2111 IN SIGMAR SUBROUTINE) 
CALL LEAVE 
STOP 

13 RETURN 
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SUBROUTINE CHISQ 

IF DIVIDE CHECK 10,11 

10 WRITE OUTPUT TAPE 6 , 1 1 

110 PORMAT(59H DIVIDE CHECK TRIGGER FOUND ON AT START OF CHISQ SUBROUT 
1INE) 

CALL LEAVE 
STOP 

11 ISPILL=0 
JSPILL=0 
CHI2ST=0 
CHI2PT=0 

DO 20 J = 1,JMAX 

CHI2S ( J)= ( (SGMATH( J)-SGMAEX( J ) ) /DSGMEX( J ) ) * * 2 . 
CHI2P ( J)= ( (POLTH( J) -POLEX( J ) ) /DPOLEX( J ) ) * * 2 . 
IF DIVIDE CHECK 14,15 

14 WRITE OUTPUT TAPE 6 , 1 1 4 , J 

114 FORMAT(40H DIVISOR DSGMEX OR DPOLEX IS ZERO FOR J = 13,20H IN CHISQ 
lSUBROUTDNE) 
CALL LEAVE 
STOP 

15 CHI2ST=CHI2ST + CHI2S(J) 
CHI2 ( J) = CHI25 ( J)+CHI2P ( J ) 

2 CHI2PT=CHI2PT+CHI2P ( J ) 
CHI2T=CHI2ST+CHI2PT 
IF (ISPILL) 16,17,16 

16 WRITE OUTPUT TAPE 6 , 1 1 6 , ISPILL 

116 FORMAT(23H UNDERFLOW OCCURRED AT 16,20H IN CHISQ SUBROUTINE) 

17 IF(ISPILL) 18,19,18 

18 WRITE OUTPUT TAPE 6,1 18, JSPILL 

118 FORMAT(22H OVERFLOW OCCURRED AT 16,20H IN CHISQ SUBROUTINE) 
CALL LEAVE 
STOP 
19 RETURN 
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SUBROUTINE 0UTPT4 

NPGS=0 

CALL SKIP(K,NPGS,NUMRUN) 
WRITE OUTPUT TAPE 6 ,245 .NUMPRG 
245 FORMAT ( 1 6HDPROGRAM NUMBER 15) 

DO 8 1 = 1,13 

WRITE OUTPUT TAPE 6 ,250 , 1 , (KTRL( I ) ) 
250 FORMAT (6H KTRL(I2 ,2H) = I2) 
8 CONTINUE 

WRITE OUTPUT TAPE 6,12 
1 2 FORMAT ( 1 1 HOB ASIC DATA) 
FKAYA=FKAY*A 
FKAYB=FKAY*BG 

WRITE OUTPUT TAPE 6 , 14 ,FMI,FMB,ELAB, ZZ , V,W, A,RO, VS ,WS,RC,BG,RG, 
14 FORMAT(7H0MSUBI=E16.9 ,10H MSUBB=E16 .9 , 10H ELAB=E16 .9 , 10H 

1 ZZP=E16.9/7H0 V=E16.9,10H V\=E16.9,10H A=E16.9, 

210H RO=E16.9/7H0 VS=E16.9,10H WS=E16.9,36H 

3 ROE16.9/5 9H0 

4 BG=E16.9,10H RG=E16.9) 

WRITE OUTPUT TAPE 6 , 16 ,RHOBN,RHOBC,RHOBNG,ECM,ETA,FKAY,FKAYA,FKAYB 
1 6 FORMAT( 7H0RHOBN=E16 . 9 , 1 OH RHOBC=E16 . 9 , 1 OH RHOBNG=E16 . 9 , 1 OH 

1 EQVfc=E16.9/7H0 ETA=E16 .9 , 1 OH K=E16.9/10H KA=E16.9, 

210H KB=E16.9) 

KT=KTRL( 7) +KTRL( 8 ) +KTRL( 9 ) +KTRL ( 1 ) 
IF (KT) 13,1818,13 
1 3 WRITE OUTPUT TAPE 6,150 ,HA,RMA, FN1A , FN2A ,PMA, HB ,RMB, FN1B , FN2B ,PMB 
150 FORMAT(7H0 HA=E16.9,7H RMA=E16.9,7H N1A=E16.9,7H N2A=E16.9 
1,7H PMA=E16.9/7H HB=E16.9,7H RMB=E16.9,7H N1B=E16.9,7H 

2N2B=E16.9 ,7H PMB=E16.9) 
1818 WRITE OUTPUT TAPE 6 , 1 8 ,RHOMAX,LMAXM 
1 8 FORMAT ( 1 7H0MTEGRATION DATA/ 8H0RHOMAX=E 1 6 . 9 , 1 OH IMAX1VH 5 ) 

WRITE OUTPUT TAPE 6 ,220 ,NMAX 
220 FORMAT (6H0NMAX=I5) 

WRITE OUTPUT TAPE 6,24 
24 FORMAT (6H0RHOIN) 
NOLINE=50 
K=20 

DO 40 1 = 1, NMAX,6 
IF(K-NOLINE) 30,29,29 

29 CALL SKIP(K,NPGS,NUMRUN) 

30 AT=XMNOF(I+5,NMAX) 
K=K+1 

WRITE OUTPUT TAPE 6,32, (RHOM( J ) , J=I ,M) 
32 FORMAT(lH E19.9 ,5E20.9) 

40 CONTINUE 

WRITE OUTPUT TAPE 6 ,41 

41 FORMAT (7H0DRHOIN) 
DO 60 I = 1,NMAX,6 
IF(K-NOLINE) 45,43,43 

43 CALL SKIP(K,NPGS,NUMRUN) 
45 AT=XATCMOF(I+5,NMAX-l) 

K=K+1 

WRTTE OUTPUT TAPE 6,32, (DRHOIN( J ) , J=I ,M) 
60 CONTINUE 
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WRITE OUTPUT TAPE 6,118 ,SGMRTH 
118 PORMAT(12H0SIGMAR(TH)=E16.9) 
15 IF(KTRL(2)-1) 1900,20,1900 

2 WRITE OUTPUT TAPE 6,119, CHI2ST , CHI2PT , CHI2T 
119 FORMAT (25HOSUM OF CHI SQUARE SIGMA=E16.9/2 3H0SUM OF CHI SQUARE PO 
1L=E16.9/2 5H0SUM OF CHI SQUARE TOTAL=E16.9) 
21 CALL SKIP(K,NPGS,NUMRUN) 

WRITE OUTPUT TAPE 6,200 
200 FORMAT (113H THETA SIGMATH SIG-SIGC 

1 POL TH SIGMA EX POL EX) 

DO 90 I = 1,JMAX 
IF (K-NOLINE) 75,70,70 
70 CALL SKIP(K,NPGS,NUMRUN) 
75 K=K+1 

WRITE OUTPUT TAPE 6,32 ,THETAD( I ) ,SGMATH( I ) ,SRATIO( I ) ,POLTH( I ) , 
1SGMAEX( I ) ,POLEX( I ) 
90 CONTINUE 
GOTO 299 
1900 CALL SKIP (K,NPGS,NUMRUN) 

WRITE OUTPUT TAPE 6 ,1905 
1905 FORMAT (120H THETA SIGMATH 

1 SIG-SIGC POL TH 

2) 

DO 1920 I = 1,JMAX 

IF (K-NOLTNE) 1910,1908,1908 
1908 CALL SKIP (K,NPGS,NUMRUN) 
1910 K=K+1 

WRITE OUTPUT TAPE 6 , 1919 ,THETAD( I ) ,SGMATH( I ) ,SRATIO(I) ,POLTH(I) 

1919 FORMAT (1H E20 . 9 , 3E30 . 9) 

1920 CONTINUE 

299 IF(KTRL(6)-1) 300,121,300 

300 IF(KTRL(12)-1) 25,1700,25 

1700 CALL SKIP(K,NPGS,NUMRUN) 
WRITE OUTPUT TAPE 6,1701 

1701 FORMAT (92H RHO( I ) FFCR FFCI 
1 FFSR FFSI) 

DO 1709 I = 1,ILAST 

IF (K-NOLTNE) 1703,1702,1702 

1702 CALL SKIP (K,NPGS,NUMRUN) 

1703 WRITE OUTPUT TAPE 6,158,RHO(I) ,FFCR(I) ,FFCI(I) ,FFSR(I) ,FFSI(I) 
158 P0RMAT(1H 5E20.9) 

1709 CONTINUE 
25 IF(KTRL(2)-1) 23,22,23 
22 CALL SKIP (K,NPGS,NUMRUN) 
WRITE OUTPUT TAPE 6,95 

95 FORMAT(120H THETA DSIGMA EX DPOL EX 

1 CHI SQUARE SIGMA CHI SQUARE POL CHI SQUARE TOTAL ) 

DO 120 J = 1,JMAX 
IF(K-NOLrNE) 97,96,96 

96 CALL SKIP(K,NPGS,NUMRUN) 

97 K=K+1 

WRITE OUTPUT TAPE 6,32 ,THETAD( J ) ,DSGMEX( J ) ,DPOLEX( J ) * CHI2S ( J ) , 
1CHI2P(J))CHI2(J) 
120 CONTINUE 
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23 CALL SKIP(K,NPGS,NUMRUN) 
1623 WRITE OUTPUT TAPE 6,1150 

1150 FORMAT (120H L REALC(L+l/2) IMA 

1G C(L + l/2) REAL C(L-l/2) DvlAG C(L-l/2) 

2) 

DO 160 L=1,LMAX 
IF (K-NOLrNE) 155,153,153 
153 CALL SKIP (K,NPGS,NUMRUN) 
155 K=K+1 
L1=L-1 

WRITE OUTPUT TAPE 6 , 1 156 ,L1 ,CR1(L) ,CI1(L) ,CR2(L) ,CI2(L) 
1156 FORMAT (1H 111 ,E30.9 ,3E25.9) 
160 CONTINUE 
121 RETURN 
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SUBROUTINE SKIP(K,NPGS,NUMRUN) 

NPGS=^PGS+1 

WRITE OUTPUT TAPE 6,1510, (NUMRUN( I ) , I = 1 ,5) ,NPGS 
1510 KKMAT(12HlRUNNUVEHtI2 ,1H-I2 ,1H-I4 ,3H -13, 3H -I3,79H 

1 PA 

2GE 15/) 

K=0 

RETURN 



SUBROUTINE LEAVE 
CALL PDUMP(A,ZZ) 
CALL CTRL4 
RETURN 



91 





CARDS 


COLUMN 




FAP 






COUNT 


43 


■SPILL 


SUBROUTINE 




ENTRY 


SPILL 


SPILL 


STZ* 


1,4 




STZ* 


2,4 




STZ 







CAL 


1,4 




STA 


AA41 




STA 


AA36 




CAL 


2,4 




STA 


AA31 




CLA* 


3,4 




STO 


AA45 




CLA* 


4,4 




STO 


AA46 




CAL 


AA47 




SLW 


8 




TRA 


5,4 


AA16 


LDI 







LFT 


4 




TRA 


AA36 




LFT 


16 




TRA 


AA24 




TRA* 





AA24 


LNT 


1 




TRA* 







CAL 







SUB 


AA35 




LLS 


18 


AA31 


STD 


AA31 




CLA 


AA46 




LDQ 


AA46 




TRA* 





AA35 


HTR 


1 


AA36 


CLA 


AA36 




TNZ 


AA42 




CAL 







SUB 


AA35 




LLS 


18 


AA41 


STD 


AA41 


AA42 


CLA 


AA45 




LDQ 


AA45 




TRA* 





AA45 


HTR 





AA46 


HTR 





AA47 


TRA 
END 


AA16 



STORE ZERO IN JSPILL 
STORE ZERO IN ISPILL 
STORE ZERO IN LOCATION 00000 

SET ADDRESS AA41 , 
AA36 TO JSPILL 
SET ADDRESS AA31 
TO ISPILL 

SET GOMVONr STORAGE 

SET OOVIVCN STORAGE 

PLACE TRANSFER 
INSTRUCTION IN LOCATION 8 

EXIT TO MAIN PROGRAM 
ENTRY IN CASE OF OVER-OR UNDERFLOW 
TEST FOR OVERFLOW 
TRANSFER IN CASE OF OVERFLOW 

TRANSFER IN CASE OF UNDERFLOW 
TRANSFER TO MADN PROGRAM NO UFLOW 
TEST FOR UNDERFLOW 
UNDERFLOW IN AC ONLY 
PLACE LOCATION AT WHICH 
UNDERFLOW OCCURRED IN AC 
SHIFT LEFT 18 
STORE IN ISPILL 
SET AC, MQ WITH 

SPECIFIED CONSTANTS 
EXIT TO MALN PROGRAM 
CONSTANT 

TEST IF JSPILL ZERO 
TRANSFER IN CASE JSPILL NON-ZERO 
PLACE LOCATION AT WHICH OVERFLOW OCCURRED 
IN AC 
SHIFT LEFT 18 
STORE IN JSPILL 
SET AC,MQ WITH SPECIFIED CONSTANTS 

EXIT TO MAIN PROGRAM 
GOMVCN STORAGE 
GOMVON STORAGE 
INSTRUCTION TO BE INSERTED AT LOC . 8 



VII. Typical Input and Output 



A. Input Data for Protons against Copper at 9.75 MeV 



3 +0.62500000 -01 +0.40750000 +01 +0.00000000 +00 
22 +0.25000000 +00 +0.00000000 +00 -0.16000000 +00 

1960 10 +0.33390000 +01 -0.20000000 +00 

32 +0.00000000 +00 +0.00000000 +00 

+0.15200000 +02 +0.33560000 +01 -0.17000000 +00 

4 +0.20300000 +02 +0.37570000 +01 -0.17000000 +00 

+0.25400000 +02 +0.38570000 +01 +0.00000000 +00 

1 +0.28000000 +02 +0.00000000 +00 -0.10000000 +00 
1 +0.30400000 +02 +0.38460000 +01 +0.00000000 +00 

+0.33000000 +02 +0.00000000 +00 +0.10000000 -01 

1 +0.35500000 +02 +0.37570000 +01 +0.00000000 +00 
+0.39000000 +02 +0.00000000 +00 +0.20000000 +00 
+0.40600000 +02 +0.39800000 +03 +0.00000000 +00 
+0.43000000 +02 +0.35500000 +02 +0.00000000 +00 
+0.45600000 +02 +0.16700000 +02 +0.00000000 +00 
+0.47000000 +02 +0.10000000 +30 +0.13000000 +00 
+0.507000000+02 +0.90800000 +01 +0.00000000 +00 

+0.51500000 +02 +0.10000000 +30 +0.70000000 -01 

1 +0.54000000 +02 +0.53800000 +01 +0.00000000 +00 
fO. 10000000 +01 +0.55700000 +02 +0.10000000 +30 -0.20000000 -01 
fO. 64000000 +02 +0.57000000 +02 +0.37300000 +01 +0.10000000 +30 
fO. 97500000 +01 +0.60000000 +02 +0.10000000 +30 +0.10000000 +30 
fO. 29000000 +02 +0.60800000 +02 +0.19100000 +01 +0.10000000 +30 
fO. 12000000 +01 +0.65500000 +02 +0.10000000 +30 +0.30000000 -01 
fO. 62000000 +02 +0.65800000 +02 +0.91500000 +00 +0.10000000 +30 
fO. 85000000 +01 +0.69000000 +02 +0.10000000 +30 +0.40000000 -01 
fO. 12000000 +01 +0.70800000 +02 +0.10000000 +30 +0.10000000 +30 
fO. 52000000 +00 +0.75500000 +02 +0.49600000 +00 +0.30000000 -01 
-0.40000000 +01 +0.75900000 +02 +0.10000000 +30 +0.10000000 +30 
fO. 00000000 +00 +0.80900000 +02 +0.10000000 +30 +0.30000000 -01 
fO. 00000000 +00 +0.85900000 +02 +0.25800000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.86000000 +02 +0.10000000 +30 +0.30000000 -01 
fO. 00000000 +00 +0.90900000 +02 +0.16300000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.95500000 +02 +0.10000000 +30 +0.40000000 -01 
fO. 00000000 +00 +0.95900000 +02 +0.13400000 +00 +0.40000000 -01 
fO. 00000000 +00 +0.10000000 +03 +0.10000000 +30 +0.10000000 +30 
fO. 00000000 +00 +0.38650000 +04 +0.13400000 +00 +0.40000000 -01 
fO. 00000000 +00 +0.97340000 +03 +0.15000000 +00 +0.30000000 -01 
fO. 00000000 +00 +0.42470000 +03 +0.15400000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.00000000 +00 +0.10000000 +30 +0.50000000 -01 
fO. 00000000 +00 +0.22690000 +03 +0.15400000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.00000000 +00 +0.10000000 +30 +0.40000000 -01 
fO. 00000000 +00 +0.13460000 +03 +0.15000000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.00000000 +00 +0.10000000 +30 +0.60000000 -01 
fO. 00000000 +00 +0.82920000 +02 +0.00000000 +00 +0.10000000 +30 
fO. 00000000 +00 +0.00000000 +00 +0.00000000 +00 +0.10000000 +30 

1 +0.47660000 +02 +0.00000000 +00 +0.10000000 +30 

1 +0.00000000 +00 -0.20000000 -01 +0.60000000 -01 

1 +0.22870000 +02 +0.00000000 +00 +0.10000000 +30 

1 +0.00000000 +00 +0.10000000 -01 +0.50000000 -01 

1 +0.00000000 +00 +0.00000000 +00 +0.10000000 +30 

1 +0.12410000 +02 -0.30000000 -01 +0.60000000 -01 

3 +0.00000000 +00 +0.00000000 +00 100 

fO. 62500000 -01 +0.00000000 +00 -0.60000000 -01 

fO. 50000000 +00 +0.64560000 +01 +0.00000000 +00 

fO. 10000000 +02 +0.00000000 +00 -0.10000000 +00 
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B. Output Listing 

KUNNUVOafc 2-40-1961 - 1 
PAGE 1 

PROGRAM NUMBER 4 

KTRL( 1) = 
KTRL( 2)=1 
KTRL( 3) = 1 
KTRL( 4)=0 
KTRL( 5) = 1 
KTRL( 6)=0 
KTRL( 7) = 
KTRL( 8)=0 
KTRL( 9)=0 
KTRL(10) = 
KTRL(11) = 
KTRL(12) = 
KTRL(13) = 1 



MSUBB= 0.639999993E 02 
Vf= 0.849999994E 01 
WS= 0. 

RHOBG= 0.323980615E 01 
K= 0.674959674E 00 



IMAXIVfc 



10 



BASIC DATA 

MSUB1= 0.099999994E 01 

V= 0.619999997E 02 
VS=-0.399999999E 01 
9 
RHOBN= 0.393980615E 01 
ETA= 0.146788672E 01 
INTEGRATION DATA 
RHOMAJfc 0.099999994E 02 
NMAJfc 3 

RHOIN 

0.625000000E-01 
DRHOIN 

0.625000000E-01 0.250000000E-00 

SIGMAR(TH)= 0.668857820E 02 
SUM OF CHI SQUARE SIGMA= 0.587550342E 02 
SUM OF CHI SQUARE POL= 0.999665476E 02 
SUM OF CHI SQUARE TOTAU= 0.158721581E 03 



ELAB= 0.974999994E 01 
A= 0.519999996E 00 

BG= 0. 
RHOBNG= 0. 

KA= 0.350979023E-00 



ZZP= 0.289999999E 02 

RO= 0.119999997E 01 

RC= 0.119999997E 01 

RG= 0. 

ECM= 0.959999986E 01 

KB= 0. 



0.500000000E 00 



0.099999994E 02 
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VIII. Further Subroutines and Programs in 

Preparation 



The following subroutines are presently being prepared at UCLA: 

Subroutine TV 

This subroutine is designed to output on CRT and on film various required curves 
such as (t(6) vs 6, a($)/a c ($) vs 6, P(6) vs 6. 

Subroutine RHOBEG 

This subroutine will make use of the quantities IIN(L) to allow the numerical integra- 
tions to start at different values of p depending upon £ in order to speed up the numerical 
integration. 

Subroutine FLUX 

This subroutine will if desired compute the normalized total wave functions, the scat- 
tered flux j, the divergence and the curl of j at specified values of p and 6. 

All the above subroutines will of course require some modification of the basic program. 
The following programs are presently being prepared at UCLA: 

Program SCAT 3 

This program will be similar to program SCAT 4 except that it will treat incident and 
target particles of zero spin, thus speeding up the calculation for that case. 

Program SCAT 5 

This is a modified version of program SCAT 4 offering a simplified input and using 
only as many £'s as may be significant in the C/s calculations. 

Program SCAT K 

This is a modified version of program SCAT 4 designed to analyze the scattering of 
K-mesons against complex nuclei, including the use of an approximate Klein-Gordon 
equation, relativistic kinematic corrections, and averaging of the cross sections over angles, 
energies, and representative nuclei. 

Program SCAT 6 

This is a modified version of program SCAT 4 designed to calculate cross sections and 
polarization of spin 1 particles scattered by spin targets. 

Program SEEK 4 

This is a program designed to search automatically the parameter space so as to 

• • • 2 
minimize \ ■ 
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